BSTS Selected Models

rm(list=ls());gc()

##          used (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells 464418 24.9    1001458  53.5   650900  34.8
## Vcells 860757  6.6   17698584 135.1 20029491 152.9

unlink('Consolidacion_BDs_cache', recursive = TRUE)
unlink('Consolidacion_BDs_FINAL_cache', recursive = TRUE)
#load(url("https://drive.google.com/uc?export=download&id=1zLfLpJjGCfMSfwrtTnaOdBSgGg2ZCqjR"))
load(paste0(getwd(),"/","Procesos hasta 4_2.RData"))

#xaringan::inf_mr()

if(isTRUE(getOption('knitr.in.progress'))==T){
    clus_iter=30000
} else {
  input <- readline('¿Are you gonna run the dataset with the whole iterations? (Si/No): ')
  if(input=="Si"){
    clus_iter=30000
  } else {
    clus_iter=1000
  }
}

#arriba puse algunas opciones para que por defecto escondiera el código
#también cargue algunos estilo .css para que el texto me apareciera justificado, entre otras cosas.
local({r <- getOption("repos")
       r["CRAN"] <- "http://cran.r-project.org" 
       options(repos=r)
})

`%>%` <- magrittr::`%>%`
copy_names <- function(x,row.names=FALSE,col.names=TRUE,dec=",",...) {
  if(class(ungroup(x))[1]=="tbl_df"){
        if(options()$OutDec=="."){
            options(OutDec = dec)
            write.table(format(data.frame(x)),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ".")
          return(x)
        } else {
            options(OutDec = ",")
            write.table(format(data.frame(x)),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ",")
          return(x)    
        }
  } else {
        if(options()$OutDec=="."){
            options(OutDec = dec)
            write.table(format(x),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ".")
          return(x)
        } else {
            options(OutDec = ",")
            write.table(format(x),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ",")
          return(x)       
  }
 }
}  

unlink('Consolidacion_BDs_cache', recursive = TRUE)
if(!require(pacman)){install.packages("pacman")}

pacman::p_unlock(lib.loc = .libPaths()) #para no tener problemas reinstalando paquetes
knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE
)
#dejo los paquetes estadísticos que voy a utilizar

if(!require(plotly)){install.packages("plotly")}
if(!require(lubridate)){install.packages("lubridate")}
if(!require(htmlwidgets)){install.packages("htmlwidgets")}
if(!require(tidyverse)){install.packages("tidyverse")}
if(!require(gganimate)){install.packages("gganimate")}
if(!require(readr)){install.packages("readr")}
if(!require(stringr)){install.packages("stringr")}
if(!require(data.table)){install.packages("data.table")}
if(!require(DT)){install.packages("DT")}
if(!require(ggplot2)){install.packages("ggplot2")}
if(!require(lattice)){install.packages("lattice")}
if(!require(forecast)){install.packages("forecast")}
if(!require(zoo)){install.packages("zoo")}
if(!require(panelView)){install.packages("panelView")}
if(!require(janitor)){install.packages("janitor")}
if(!require(rjson)){install.packages("rjson")}
if(!require(estimatr)){install.packages("estimatr")} 
if(!require(CausalImpact)){install.packages("CausalImpact")}
if(!require(textreg)){install.packages("textreg")}
if(!require(sjPlot)){install.packages("sjPlot")}
if(!require(foreign)){install.packages("foreign")}
if(!require(tsModel)){install.packages("tsModel")}
if(!require(lmtest)){install.packages("lmtest")}
if(!require(Epi)){install.packages("Epi")}
if(!require(splines)){install.packages("splines")}
if(!require(vcd)){install.packages("vcd")}
if(!require(astsa)){install.packages("astsa")}
if(!require(forecast)){install.packages("forecast")}
if(!require(MASS)){install.packages("MASS")}
if(!require(ggsci)){install.packages("ggsci")}
if(!require(Hmisc)){install.packages("Hmisc")}
if(!require(compareGroups)){install.packages("compareGroups")}
if(!require(dplyr)){install.packages("dplyr")}
if(!require(ggforce)){install.packages("ggforce")}
if(!require(imputeTS)){install.packages("imputeTS")}
if(!require(doParallel)){install.packages("doParallel")}
if(!require(SCtools)){install.packages("SCtools")}
if(!require(MSCMT)){install.packages("MSCMT")}
if(!require(ISOweek)){install.packages("ISOweek")}
if(!require(gridExtra)){install.packages("gridExtra")}
if(!require(cowplot)){install.packages("cowplot")}
if(!require(grid)){install.packages("grid")}
if(!require(ArCo)){install.packages("ArCo")}

# Calculate the number of cores
no_cores <- detectCores() - 1
cl<-makeCluster(no_cores)
registerDoParallel(cl)

Sys.setlocale(category = "LC_ALL", locale = "english")

## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:PARA GENERAR MATRICES DE PREDICCION BSTS#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

GetPosteriorStateSamples <- function(bsts.model) {
  # Returns a matrix of simulated values from the marginal posterior
  # distribution for the sum of all state variables.
  #
  # Args:
  #   bsts.model: A fitted model returned by \code{bsts()}.
  #
  # Returns:
  #   matrix [number of post-burn-in MCMC samples] x [time points]
  
  # Get state contributions (e.g., 1000 samples x 2 states x 365 time pts),
  # discarding burn-in samples (=> 900 x 2 x 365)
  burn <- SuggestBurn(0.1, bsts.model)
  #assert_that(burn > 0)
  state.contributions <- bsts.model$state.contributions[-seq_len(burn), , ,
                                                        drop = FALSE]
  
  # Sum across states, call it 'state.samples' (=> 900 x 365)
  state.samples <- rowSums(aperm(state.contributions, c(1, 3, 2)), dims = 2)
  return(state.samples)
}
ComputeResponseTrajectories <- function(bsts.model) {
  # Generates trajectories of the response variable. A trajectory is a simulated
  # time series drawn from the posterior predictive distribution over the data.
  # This function differs from GetPosteriorStateSamples(). The latter returns
  # the posterior mean of the response. This function returns the actual value
  # (posterior mean + observation noise).
  #
  # Args:
  #   bsts.model: A model object as returned by \code{bsts()}.
  #
  # Returns:
  #   matrix [number of post-burn-in MCMC samples] x [time points]
  
  # Get posterior state samples
  state.samples <- GetPosteriorStateSamples(bsts.model)
  
  # Get observation noise standard deviation samples
  burn <- SuggestBurn(0.1, bsts.model)
  #assert_that(burn > 0)
  sigma.obs <- bsts.model$sigma.obs[-seq_len(burn)]  # e.g., 900
  
  # Sample from the posterior predictive density over data
  n.samples <- dim(state.samples)[1]  # e.g., 900 x 365
  obs.noise.samples <- matrix(rnorm(prod(dim(state.samples)), 0, sigma.obs),
                              nrow = n.samples)
  y.samples <- state.samples + obs.noise.samples
  return(y.samples)
}


#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
aicbic_plm <- function(object, criterion) {
  
  
  # object is "plm", "panelmodel" 
  # Lets panel data has index :index = c("Country", "Time")
  
  sp = summary(object)
  
  if(class(object)[1]=="plm"){
    u.hat <- residuals(sp) # extract residuals
    df <- cbind(as.vector(u.hat), attr(u.hat, "index"))
    names(df) <- c("resid", "Country", "Time")
    c = length(levels(df$Country)) # extract country dimension 
    t = length(levels(df$Time)) # extract time dimension 
    np = length(sp$coefficients[,1]) # number of parameters
    n.N = nrow(sp$model) # number of data
    s.sq  <- log( (sum(u.hat^2)/(n.N))) # log sum of squares
    
    # effect = c("individual", "time", "twoways", "nested"),
    # model = c("within", "random", "ht", "between", "pooling", "fd")
    
    # I am making example only with some of the versions:
    
    if (sp$args$model == "within" & sp$args$effect == "individual"){
      n = c
      np = np+n+1 # update number of parameters
    }
    
    if (sp$args$model == "within" & sp$args$effect == "time"){
      T = t
      np = np+T+1 # update number of parameters
    }
    
    if (sp$args$model == "within" & sp$args$effect == "twoways"){
      n = c
      T = t
      np = np+n+T # update number of parameters
    }
    aic <- round(       2*np  +  n.N * (  log(2*pi) + s.sq  + 1 ),1)
    bic <- round(log(n.N)*np  +  n.N * (  log(2*pi) + s.sq  + 1 ),1)
    
    if(criterion=="AIC"){
      names(aic) = "AIC"
      return(aic)
    }
    if(criterion=="BIC"){
      names(bic) = "BIC"
      return(bic)
    }
  }
}

1 Dataset Compile

library(compareGroups)
#pvals <- compareGroups::getResults(table, "p.overall")
table2 <- compareGroups::compareGroups(did ~ cons_total+ cons_trauma+ cons_resp+ cons_circ+ hosp_total+ hosp_trauma+ hosp_resp+ hosp_circ+ rate+ rate_resp, #27
                       method = c(cons_total=2, cons_trauma=2, cons_resp=2, cons_circ=2,  hosp_total=2,
                                  hosp_resp=2, hosp_circ=2, hosp_trauma=2,rate=2, rate_resp=2),
                      data = data15a64_rn_ratio,
                      include.miss = T,
                      var.equal=T)
                      #,
                      #subset = age == "15-64")
#p.adjust(pvals, method = "BH")
restab <- createTable(table2, show.n = F, show.p.overall = F)
 compareGroups::export2md(restab, first.strip = T, show.all = T,hide.no = "no",header.labels = c(`0`= "Pre- Oct 18, 2020", `1`="Post Oct 18, 2020"), size=12,col.names= c("", "Previous to social
protests","During social protests"),format="html",position="center",caption= "Table 1. Median descriptive table of emergency department consultations and hospitalizations, pre and port Ocober’s 2019 social protests in Chile")%>%#,p.overall = "p-value"
  kableExtra::add_footnote(paste0("Note. Percentiles 25 and 75 in brackets."), notation = "none")%>%
  kableExtra::scroll_box(width = "100%", height = "375px")

Table 1. Median descriptive table of emergency department consultations and hospitalizations, pre and port Ocober’s 2019 social protests in Chile

	Previous to social protests	During social protests
	N=252	N=10
Total Consultations	3137 [2924;3361]	2854 [2754;2898]
Trauma Consultations	802 [728;888]	786 [752;801]
Respiratory Consultations	143 [120;183]	96.0 [77.5;103]
Circulatory Consultations	102 [87.0;125]	90.5 [87.5;95.8]
Total Hospitalizations	288 [268;311]	298 [281;332]
Trauma Hospitalizations	60.0 [52.0;67.0]	81.5 [77.5;89.8]
Respiratory Hospitalizations	19.5 [15.0;23.2]	22.0 [15.5;24.0]
Circulatory Hospitalizations	29.0 [23.0;36.0]	35.5 [30.5;42.0]
Rate of Trauma Hospitalizations per Trauma Consultations (x1,000 population)	73.0 [64.0;86.2]	102 [84.2;113]
Respiratory Hospitalizations per Respiratory Consultations (x1,000 population)	131 [107;160]	233 [189;270]
Note. Percentiles 25 and 75 in brackets.

2 Bayesian Structural Time Series Models

The model of Bayesian Structural Time Series (BSTS) with Causal Impact Analysis contains a local level component (observation equation linking observed data to a state vector) and a seasonal component (state equation that describes how the state vector evolves over time). Also it can contain other predictor series. The method uses a Markov chain Monte Carlo (MCMC) sampling of the structural time series given the observed data and the priors. We run MCMC simulations with 10,000 draws.

Previously, we selected those specifications on time series that showed less cumulative errors for each outcome. Here, we defined the final models by controlling only with circulatory system hospitalizations. We discarded other covariates, because we suspected that respiratory hospitalizations or consultations, trauma hospitalizations or consultations, and total consultations or total hospitalizations (composed by respiratory and trauma) may also be affected by social protests.

For all the models, we used the term prior.level.sd to .1 instead of the default .01, assuming an unstable dataset with high residual volatility, although may give rise to wide prediction intervals.

rm(list=ls()[startsWith(ls(),"impact")])
rm(list=ls()[startsWith(ls(),"model")])
#We did not include time-varying regression coefficients, because this may leads to overspecification when used in combination with a time-varying local trend or level, in which case a static regression is safer.

2.1 Trauma Consultations

The model that had lower cumulative errors assumed a Random-Walk model and a Studentized Distribution with control variables.

#The response variable (i.e., the first column in data) may contain missing values (NA), but covariates (all other columns in data) may not. If one of your covariates contains missing values, consider imputing (i.e., estimating) the missing values; if this is not feasible, leave the regressor out.

# Model 2
ss2d <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2d <- AddLocalLevel(ss2d, 
                      c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"cons_trauma"])), rep(NA,10))
                      ) #
# Add weekly seasonal
ss2d <- AddSeasonal(ss2d, 
                    c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"cons_trauma"])), rep(NA,10)),
                    nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2d <- AddSeasonal(ss2d, 
                    c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"cons_trauma"])), rep(NA,10)), 
                    nseasons = 12, season.duration =4) #years
#ss2 <- AddAutoAr(ss2, y = data15a64_rn_causal$log_hosp_trauma, lags = 1) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model2d1_cons_trauma <-  bsts(c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"cons_trauma"])), rep(NA,10)),
                   state.specification = ss2d,#A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter, 
              # burn = 500, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)

## =-=-=-=-= Iteration 0 Tue Feb 09 07:43:58 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 07:44:13 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 07:44:28 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 07:44:43 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 07:44:58 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 07:45:13 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 07:45:28 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 07:45:44 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 07:45:59 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 07:46:14 2021
##  =-=-=-=-=

#,
#               dynamic.regression=T)
#plot(model2d1_cons_trauma, main = "Model of Trauma Consultations (Circulatory System Hospitalizations as Control Var)")
#plot(model2d1_cons_trauma, "components")

impact3d1_cons_trauma <- CausalImpact(bsts.model = model2d1_cons_trauma,model.args = list(prior.level.sd=.1, dynamic.regression=T),
                       post.period.response = as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==1),"cons_trauma"])))
#plot(impact3d1_cons_trauma, "original") 

burn2d1 <- SuggestBurn(0.1, model2d1_cons_trauma)

#summary(impact3d1_cons_trauma)
#summary(impact3d1) 
##summary(impact3d1,"report")

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
model2d1_cons_trauma_resp_traj<-(ComputeResponseTrajectories(model2d1_cons_trauma))
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

plot(impact3d1_cons_trauma)+
  xlab("Date")+
  ylab("Consultations")+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn_ratio),12),262)),
                     labels=as.character(unlist(data15a64_rn_ratio[c(seq(1,nrow(data15a64_rn_ratio),12),262),"date"])))+
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=11),
        plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  labs(caption="Note. The first panel shows the data and a counterfactual prediction for the post-treatment period (Blue dashed line);\nThe second panel shows the difference between observed data and counterfactual predictions;\nThe third panel adds up the pointwise contributions from the second panel;\nBlue area= Prediction intervals.")

2.2 Respiratory Consultations

The model that had lower cumulative errors assumed a Random-Walk model and a Studentized Distribution with control variables.

#The response variable (i.e., the first column in data) may contain missing values (NA), but covariates (all other columns in data) may not. If one of your covariates contains missing values, consider imputing (i.e., estimating) the missing values; if this is not feasible, leave the regressor out.

# Model 2
ss2d <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2d <- AddLocalLevel(ss2d,c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"cons_resp"])), rep(NA,10))
                      ) #
# Add weekly seasonal
ss2d <- AddSeasonal(ss2d, 
                    c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"cons_resp"])), rep(NA,10)), 
                    nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2d <- AddSeasonal(ss2d, 
                    c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"cons_resp"])), rep(NA,10)), 
                    nseasons = 12, season.duration =4) #years
#ss2 <- AddAutoAr(ss2, y = data15a64_rn_causal$hosp_trauma, lags = 1) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).

model2d_cons_resp <- bsts(c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"cons_resp"])), rep(NA,10)),
               state.specification = ss2d, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter, 
              # burn = 500, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               model.options = BstsOptions(save.prediction.errors = TRUE),
               seed= 2125)

## =-=-=-=-= Iteration 0 Tue Feb 09 07:47:10 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 07:47:24 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 07:47:38 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 07:47:52 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 07:48:06 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 07:48:20 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 07:48:34 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 07:48:48 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 07:49:02 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 07:49:16 2021
##  =-=-=-=-=

#,
#               dynamic.regression=T)
#plot(model2d_cons_resp, main = "Model 2")
#plot(model2d_cons_resp, "components")

impact3d_cons_resp <- CausalImpact(bsts.model = model2d_cons_resp,model.args = list(prior.level.sd=.1, dynamic.regression=T),
                       post.period.response = as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==1),"cons_resp"])))
#plot(impact3d, "original") 

burn2d <- SuggestBurn(0.1, model2d_cons_resp)

#summary(impact3d_cons_resp)

2.3 Trauma Hospitalizations

The model that had lower cumulative errors assumed a Random-Walk model and a Studentized Distribution with control variables.

#The response variable (i.e., the first column in data) may contain missing values (NA), but covariates (all other columns in data) may not. If one of your covariates contains missing values, consider imputing (i.e., estimating) the missing values; if this is not feasible, leave the regressor out.

# Model 2
ss2d <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2d <- AddLocalLevel(ss2d, c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"hosp_trauma"])), rep(NA,10))) #
# Add weekly seasonal
ss2d <- AddSeasonal(ss2d, c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"hosp_trauma"])), rep(NA,10)), nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2d <- AddSeasonal(ss2d, c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"hosp_trauma"])), rep(NA,10)), 
                     nseasons = 12, season.duration =4) #years
#ss2 <- AddAutoAr(ss2, y = data15a64_rn_causal$log_hosp_trauma, lags = 1) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model2d_trauma_hosp <- bsts(c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"hosp_trauma"])), rep(NA,10)),
                state.specification = ss2d,             
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter, 
              # burn = 500, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)

## =-=-=-=-= Iteration 0 Tue Feb 09 07:50:05 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 07:50:19 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 07:50:33 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 07:50:47 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 07:51:02 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 07:51:16 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 07:51:30 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 07:51:44 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 07:51:59 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 07:52:13 2021
##  =-=-=-=-=

#,
#               dynamic.regression=T)
#plot(model2d_trauma_hosp, main = "Model 2")
#plot(model2d_trauma_hosp, "components")

impact3d_hosp_trauma <- CausalImpact(bsts.model = model2d_trauma_hosp,model.args = list(prior.level.sd=.1, dynamic.regression=T), 
                       post.period.response =as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==1),"hosp_trauma"])),
                        )
#plot(impact3d_hosp_trauma, "original") 

burn2d <- SuggestBurn(0.1, model2d_trauma_hosp)

#plot(impact3d)
#summary(impact3d_hosp_trauma) 
##summary(impact3d,"report")
#dynamic.regression Whether to include time-varying regression coefficients. In combination with a time-varying local trend or even a time-varying local level, this often leads to overspecification, in which case a static regression is safer. Defaults to FALSE.
#
#Prior standard deviation of the Gaussian random walk of the local level. Expressed in terms of data standard deviations. Defaults to 0.01, a typical choice for well-behaved and stable datasets with low residual volatility after regressing out known predictors (e.g., web searches or sales in high quantities). When in doubt, a safer option is to use 0.1, as validated on synthetic data, although this may sometimes give rise to unrealistically wide prediction intervals.
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
model2d_trauma_hosp_resp_traj<-(ComputeResponseTrajectories(model2d_trauma_hosp))
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

plot(impact3d_hosp_trauma)+
  xlab("Date")+
  ylab("Hospitalizations")+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn_ratio),12),262)),
                     labels=as.character(unlist(data15a64_rn_ratio[c(seq(1,nrow(data15a64_rn_ratio),12),262),"date"])))+
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=11),
        plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  labs(caption="Note. The first panel shows the data and a counterfactual prediction for the post-treatment period (Blue dashed line);\nThe second panel shows the difference between observed data and counterfactual predictions;\nThe third panel adds up the pointwise contributions from the second panel;\nBlue area= Prediction intervals.")

2.4 Respiratory Hospitalizations

The model that had lower cumulative errors assumed a Random-Walk model and a Studentized Distribution, but without control variables. We decided to include control variables in order to reflect most adequately the uniqueness of the impact of the social protests in respiratory hospitalizations.

#The response variable (i.e., the first column in data) may contain missing values (NA), but covariates (all other columns in data) may not. If one of your covariates contains missing values, consider imputing (i.e., estimating) the missing values; if this is not feasible, leave the regressor out.

# Model 2
ss2d <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2d <- AddLocalLevel(ss2d, c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"hosp_resp"])), rep(NA,10))
                      ) #
# Add weekly seasonal
ss2d <- AddSeasonal(ss2d, c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"hosp_resp"])), rep(NA,10)),
                    nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2d <- AddSeasonal(ss2d, c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"hosp_resp"])), rep(NA,10)), 
                    nseasons = 12, season.duration =4) #years
#ss2 <- AddAutoAr(ss2, y = data15a64_rn_causal$hosp_trauma, lags = 1) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model2d1_hosp_resp <- 
  bsts(c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"cons_trauma"])), rep(NA,10)),
               state.specification = ss2d, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter,
              # burn = 500, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)

## =-=-=-=-= Iteration 0 Tue Feb 09 07:53:05 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 07:53:20 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 07:53:35 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 07:53:50 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 07:54:05 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 07:54:20 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 07:54:35 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 07:54:50 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 07:55:05 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 07:55:19 2021
##  =-=-=-=-=

#,
#               dynamic.regression=T)
#plot(model2d1_hosp_resp, main = "Model of Respiratory Hospitalizations (Circulatory System Hospitalizations as Control Var)")
#plot(model2d1_hosp_resp, "components")

impact3d1_hosp_resp <- CausalImpact(bsts.model = model2d1_hosp_resp,model.args = list(prior.level.sd=.1, dynamic.regression=T),
                       post.period.response = as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==1),"hosp_resp"]))
                         )
#plot(impact3d1_hosp_resp, "original") 

burn2d1 <- SuggestBurn(0.1, model2d1_hosp_resp)

#summary(impact3d1_hosp_resp) 
##summary(impact3d1,"report")

plot(impact3d1_hosp_resp)+
  xlab("Date")+
  ylab("Hospitalizations")+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn_ratio),12),262)),
                     labels=as.character(unlist(data15a64_rn_ratio[c(seq(1,nrow(data15a64_rn_ratio),12),262),"date"])))+
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=11),
        plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  labs(caption="Note. The first panel shows the data and a counterfactual prediction for the post-treatment period (Blue dashed line);\nThe second panel shows the difference between observed data and counterfactual predictions;\nThe third panel adds up the pointwise contributions from the second panel;\nBlue area= Prediction intervals.")

plot(impact3d_cons_resp)+
  xlab("Date")+
  ylab("Consultations")+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn_ratio),12),262)),
                     labels=as.character(unlist(data15a64_rn_ratio[c(seq(1,nrow(data15a64_rn_ratio),12),262),"date"])))+
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=11),
        plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  labs(caption="Note. The first panel shows the data and a counterfactual prediction for the post-treatment period (Blue dashed line);\nThe second panel shows the difference between observed data and counterfactual predictions;\nThe third panel adds up the pointwise contributions from the second panel;\nBlue area= Prediction intervals.")

2.5 Trauma Hospitalizations per 1,000 consultations

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
<br>

We generated a plot of the series of the rate of hospitalizations per trauma consultations per 1,000 consultations in each year.

<br>
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
ratio_plot<-
data15a64_rn_ratio %>% 
  dplyr::select(1:4,isoweek,rate,rate_resp,rate_circ) %>% 
  pivot_longer(cols=c("rate","rate_resp","rate_circ"), names_to = "variable", values_to = "rate") %>% 
  dplyr::arrange(year_week) %>% 
  dplyr::mutate(variable=dplyr::case_when(variable=="rate"~"Trauma",
                                          variable=="rate_resp"~"Respiratory",
                                          variable=="rate_circ"~"Circulatory")) %>% 
  dplyr::mutate(label_text=paste0("<br>Date: ",year_week,"<br>Rate: ",round(rate,0),"<br>Outcome: ",variable)) %>% 
   ggplot() + #median
  facet_wrap(.~year, ncol = 1,strip.position="right") + 
  geom_line(aes(x=isoweek, y =rate, group=variable, color=variable, label= label_text)) + 
  #geom_point(aes(x=fech_week, y =count, group=variable, color=variable,shape=variable),size=1) + 
  #geom_jitter(aes(y = hosp_trauma,x = isoweek, color = year),width = 0.5, alpha = 0.3)+
  guides(color=F)+
  theme_bw() + 
  theme(strip.background  = element_blank(),
        strip.text = element_text(face="bold", size=7))+
  ylab("") + 
  theme(strip.text.x = element_text(size = 8, face = "bold"),
        legend.position = "bottom",
        plot.caption=element_text(hjust = 0),
        strip.background  = element_blank())+
  xlab("Week")+
  scale_x_continuous(
    breaks = seq(from = 1, to = 52, by =4),
    labels = seq(from = 1, to = 52, by =4),#,
    #  label = c("two", "four", "six")
  )
#+
#  scale_y_continuous(limits=c(0,.25),labels = scales::percent)
ggplotly(ratio_plot, tooltip="label_text") %>% 
  layout(showlegend = F)

The model chosen assumed a Random-Walk model and a Studentized Distribution with control variables. We selected this model because both outcomes showed less cumulative errors assuming this structure.

#https://github.com/google/CausalImpact/tree/master/R

# Model 2
ss2d <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2d <- AddLocalLevel(ss2d,c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"rate"])), rep(NA,10))
                      ) #
# Add weekly seasonal
ss2d <- AddSeasonal(ss2d, 
                    c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"rate"])), rep(NA,10)), 
                    nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2d <- AddSeasonal(ss2d, 
                    c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"rate"])), rep(NA,10)), 
                    nseasons = 12, season.duration =4) #years
#ss2 <- AddAutoAr(ss2, y = data15a64_rn_causal$hosp_trauma, lags = 1) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).

model2d_ratio <- bsts(c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"rate"])), rep(NA,10)) ~
                   as.numeric(unlist(data15a64_rn_ratio[,"rate_circ"]))                   ,
               state.specification = ss2d, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter, 
              # burn = 500, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)

## =-=-=-=-= Iteration 0 Tue Feb 09 07:56:10 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 07:56:25 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 07:56:40 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 07:56:55 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 07:57:10 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 07:57:25 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 07:57:40 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 07:57:55 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 07:58:11 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 07:58:26 2021
##  =-=-=-=-=

#,
#              
#plot(model2d_cons_resp, main = "Model 2")
#plot(model2d_cons_resp, "components")

impact3d_ratio <- CausalImpact(bsts.model = model2d_ratio,model.args = list(prior.level.sd=.1, dynamic.regression=T),
              post.period.response = as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==1),"rate"])))
#plot(impact3d, "original") 

burn2d <- SuggestBurn(0.1, model2d_ratio)

#summary(impact3d_ratio)

plot(impact3d_ratio)+
  xlab("Date")+
  ylab("Ratio")+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),12),262)),
                     labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),12),262),"date"])))+
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=11),
        plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  #scale_y_continuous(labels=scales::percent)+
  labs(caption="Note. The first panel shows the data and a counterfactual prediction for the post-treatment period (Blue dashed line);\nThe second panel shows the difference between observed data and counterfactual predictions;\nThe third panel adds up the pointwise contributions from the second panel;\nBlue area= Prediction intervals.")

<br>

We additionally looked over the matrix of MCMC draws and compared the trajectories generated in the trauma consultations and trauma hospitalizations.

<br>
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#RESPONSE TRAJECTORIES EJE- X (time), EJE-Y (mcmc.iteration)

#c(as.character(unlist(data15a64_rn[which(data15a64_rn$did==0),"date"])), rep(NA,10))
#c(as.character(unlist(data15a64_rn[,"date"])))
model2d_cons_trauma_resp_traj<-ComputeResponseTrajectories(model2d1_cons_trauma)
model2d_trauma_hosp_resp_traj<-ComputeResponseTrajectories(model2d_trauma_hosp)
model2d_prop_trauma_resp_traj<- (model2d_trauma_hosp_resp_traj/model2d_cons_trauma_resp_traj)*1000

model2d_cons_trauma_resp_traj_df<-data.frame(model2d_cons_trauma_resp_traj)
model2d_trauma_hosp_resp_traj_df<-data.frame(model2d_trauma_hosp_resp_traj)
model2d_prop_trauma_resp_traj_df<-data.frame(model2d_prop_trauma_resp_traj)

colnames(model2d_cons_trauma_resp_traj_df) <- as.character(unlist(data15a64_rn[,"date"]))
colnames(model2d_trauma_hosp_resp_traj_df) <- as.character(unlist(data15a64_rn[,"date"]))
colnames(model2d_prop_trauma_resp_traj_df) <- as.character(unlist(data15a64_rn[,"date"]))

model2d_trauma_hosp_resp_traj_df_melt<-
  model2d_trauma_hosp_resp_traj_df %>% 
      dplyr::mutate(mcmc_sample=row_number()) %>% 
      dplyr::select(mcmc_sample,everything()) %>% 
      melt(id.vars ="mcmc_sample") %>% 
      rename("date"="variable")
model2d_trauma_cons_resp_resp_traj_df_melt<-
  model2d_cons_trauma_resp_traj_df %>% 
      dplyr::mutate(mcmc_sample=row_number()) %>% 
      dplyr::select(mcmc_sample,everything()) %>% 
      melt(id.vars ="mcmc_sample") %>% 
      rename("date"="variable")
model2d_prop_trauma_resp_traj_df_melt<-
  model2d_prop_trauma_resp_traj_df %>% 
      dplyr::mutate(mcmc_sample=row_number()) %>% 
      dplyr::select(mcmc_sample,everything()) %>% 
      melt(id.vars ="mcmc_sample") %>% 
      rename("date"="variable")

model2d_rate_trauma_hosp_cons_resp_traj_df_melt<-
  model2d_trauma_hosp_resp_traj_df_melt %>% 
  dplyr::rename("trauma_hosp"="value") %>% 
  dplyr::left_join(model2d_trauma_cons_resp_resp_traj_df_melt,by = c("mcmc_sample", "date")) %>% 
  dplyr::rename("trauma_cons"="value") %>% 
  dplyr::mutate(rate_trauma=(trauma_hosp/trauma_cons)*1000)

model2d_rate_trauma_resp_traj_df_melt_final<-
model2d_rate_trauma_hosp_cons_resp_traj_df_melt %>% 
     dplyr::mutate(date=as.Date(date)) %>% 
     # dplyr::filter(date>="2019-10-21") %>% 
     dplyr::group_by(date) %>% 
    dplyr::summarise(mean = mean(rate_trauma),
                      q025 = quantile(rate_trauma, .025),
                      q975 = quantile(rate_trauma, .975),
                      n=n()) %>%
     ungroup()
#model2d_prop_trauma_resp_traj_df_melt_post<-
#    model2d_prop_trauma_hosp_cons_resp_traj_df_melt %>% 
#      dplyr::mutate(date=as.Date(date)) %>% 
#     # dplyr::filter(date>="2019-10-21") %>% 
#  dplyr::group_by(date) %>% 
#  dplyr::summarise(mean = mean(value, na.rm = TRUE),
#            sd = sd(value, na.rm = TRUE),
#            n = n()) %>%
#  dplyr::mutate(se = sd / sqrt(n),
#         lo_ci = mean - qt(1 - (0.05 / 2), n - 1) * se,
#         up_ci = mean +  (1 - (0.05 / 2), n - 1) * se) %>% 
#  ungroup()

#:#:#:#::#:PLOT#:#:#:#:#:#:
plot_comp_fit_actual_ratio<-
cbind.data.frame(model2d_rate_trauma_resp_traj_df_melt_final,
                 actual=as.numeric(unlist(data15a64_rn_ratio[,"rate"]))) %>% 
    ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=1, color="darkblue") +
      geom_line(aes(y=actual), size=1, color="black")+
  #geom_ribbon(aes(x=date, y=mean, ymin = lo_ci, ymax = up_ci), fill = "grey70")+
  geom_ribbon(aes(ymin=q025,ymax=q975),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  ylab("Rate of Hospitalizations per Consultations per 1,000 population")+
  xlab("Date")+
   theme(strip.text.x = element_text(size = 8, face = "bold"),
        legend.position = "bottom",
        plot.caption=element_text(hjust = 0),
        strip.background  = element_blank(),
        axis.text.x=element_text(angle = -90, hjust = 0),
        legend.title = element_blank())
  #scale_y_continuous(labels = scales::percent)
   # scale_x_date(labels = scales::date_format("%Y-%m"),
  #             limits = c(as.Date("2015-01-01"),as.Date("2019-12-31")),
  #               breaks = "3 months")
  #facet_zoom(xlim = c(253, nrow(actual_vs_fitted_and_forecast_arima)))
plot_comp_fit_actual_ratio+
    ggforce::facet_zoom(x = date>=as.Date("2019-10-01") & date<=as.Date("2019-12-31"),shrink=F,zoom.size = .5)+
    geom_vline(xintercept = as.Date("2019-10-21"), col = 2, lty = 2, size=1)+
  labs(caption="Note. Vertical Line, Social Protests;\nBlack Line= Actual Trend; Blue Line= Average Estimations")

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

2.6 Respiratory Hospitalizations per 1,000 consultations

The model chosen assumed a Random-Walk model and a Studentized Distribution with control variables. We selected this model because both outcomes showed less cumulative errors assuming this structure.

#https://github.com/google/CausalImpact/tree/master/R

# Model 2
ss2d <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2d <- AddLocalLevel(ss2d,c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"rate_resp"])), rep(NA,10))
                      ) #
# Add weekly seasonal
ss2d <- AddSeasonal(ss2d, 
                    c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"rate_resp"])), rep(NA,10)), 
                    nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2d <- AddSeasonal(ss2d, 
                    c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"rate_resp"])), rep(NA,10)), 
                    nseasons = 12, season.duration =4) #years
#ss2 <- AddAutoAr(ss2, y = data15a64_rn_causal$hosp_trauma, lags = 1) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).

model2d_ratio_resp <- bsts(c(as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==0),"rate_resp"])), rep(NA,10)) ~
                   as.numeric(unlist(data15a64_rn_ratio[,"rate_circ"]))                   ,
               state.specification = ss2d, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter, 
              # burn = 500, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)

## =-=-=-=-= Iteration 0 Tue Feb 09 07:59:17 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 07:59:32 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 07:59:46 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 08:00:01 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 08:00:16 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 08:00:31 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 08:00:47 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 08:01:02 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 08:01:18 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 08:01:33 2021
##  =-=-=-=-=

#,
#              
#plot(model2d_cons_resp, main = "Model 2")
#plot(model2d_cons_resp, "components")

impact3d_ratio_resp <- CausalImpact(bsts.model = model2d_ratio_resp,model.args = list(prior.level.sd=.1, dynamic.regression=T),
              post.period.response = as.numeric(unlist(data15a64_rn_ratio[which(data15a64_rn_ratio$did==1),"rate_resp"])))
#plot(impact3d, "original") 

burn2d <- SuggestBurn(0.1, model2d_ratio_resp)

#summary(model2d_ratio_resp)

plot(impact3d_ratio_resp)+
  xlab("Date")+
  ylab("Ratio")+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),12),262)),
                     labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),12),262),"date"])))+
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=11),
        plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  #scale_y_continuous(labels=scales::percent)+
  labs(caption="Note. The first panel shows the data and a counterfactual prediction for the post-treatment period (Blue dashed line);\nThe second panel shows the difference between observed data and counterfactual predictions;\nThe third panel adds up the pointwise contributions from the second panel;\nBlue area= Prediction intervals.")

<br>

We additionally looked over the matrix of MCMC draws and compared the trajectories generated in the trauma consultations and trauma hospitalizations.

<br>
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#RESPONSE TRAJECTORIES EJE- X (time), EJE-Y (mcmc.iteration)

#c(as.character(unlist(data15a64_rn[which(data15a64_rn$did==0),"date"])), rep(NA,10))
#c(as.character(unlist(data15a64_rn[,"date"])))
model2d_resp_cons_resp_traj<-ComputeResponseTrajectories(model2d_cons_resp)
model2d_resp_hosp_resp_traj<-ComputeResponseTrajectories(model2d1_hosp_resp)
model2d_prop_resp_resp_traj<- (model2d_resp_hosp_resp_traj/model2d_resp_cons_resp_traj)*1000

model2d_resp_cons_resp_traj_df<-data.frame(model2d_resp_cons_resp_traj)
model2d_resp_hosp_resp_traj_df<-data.frame(model2d_resp_hosp_resp_traj)
model2d_prop_resp_resp_traj_df<-data.frame(model2d_prop_resp_resp_traj)

colnames(model2d_resp_cons_resp_traj_df) <- as.character(unlist(data15a64_rn[,"date"]))
colnames(model2d_resp_hosp_resp_traj_df) <- as.character(unlist(data15a64_rn[,"date"]))
colnames(model2d_prop_resp_resp_traj_df) <- as.character(unlist(data15a64_rn[,"date"]))

model2d_resp_hosp_resp_traj_df_melt<-
  model2d_resp_hosp_resp_traj_df %>% 
      dplyr::mutate(mcmc_sample=row_number()) %>% 
      dplyr::select(mcmc_sample,everything()) %>% 
      melt(id.vars ="mcmc_sample") %>% 
      rename("date"="variable")
model2d_resp_cons_resp_traj_df_melt<-
  model2d_resp_cons_resp_traj_df %>% 
      dplyr::mutate(mcmc_sample=row_number()) %>% 
      dplyr::select(mcmc_sample,everything()) %>% 
      melt(id.vars ="mcmc_sample") %>% 
      rename("date"="variable")
model2d_prop_resp_resp_traj_df_melt<-
  model2d_prop_resp_resp_traj_df %>% 
      dplyr::mutate(mcmc_sample=row_number()) %>% 
      dplyr::select(mcmc_sample,everything()) %>% 
      melt(id.vars ="mcmc_sample") %>% 
      rename("date"="variable")

model2d_rate_resp_hosp_resp_traj_df_melt<-
  model2d_resp_hosp_resp_traj_df_melt %>% 
  dplyr::rename("resp_hosp"="value") %>% 
  dplyr::left_join(model2d_resp_cons_resp_traj_df_melt,by = c("mcmc_sample", "date")) %>% 
  dplyr::rename("resp_cons"="value") %>% 
  dplyr::mutate(rate_resp=(resp_hosp/resp_cons)*1000)

model2d_rate_resp_traj_df_melt_final<-
model2d_rate_resp_hosp_resp_traj_df_melt %>% 
     dplyr::mutate(date=as.Date(date)) %>% 
     # dplyr::filter(date>="2019-10-21") %>% 
     dplyr::group_by(date) %>% 
    dplyr::summarise(mean = mean(rate_resp),
                      q025 = quantile(rate_resp, .025),
                      q975 = quantile(rate_resp, .975),
                      n=n()) %>%
     ungroup()
#model2d_prop_trauma_resp_traj_df_melt_post<-
#    model2d_prop_trauma_hosp_cons_resp_traj_df_melt %>% 
#      dplyr::mutate(date=as.Date(date)) %>% 
#     # dplyr::filter(date>="2019-10-21") %>% 
#  dplyr::group_by(date) %>% 
#  dplyr::summarise(mean = mean(value, na.rm = TRUE),
#            sd = sd(value, na.rm = TRUE),
#            n = n()) %>%
#  dplyr::mutate(se = sd / sqrt(n),
#         lo_ci = mean - qt(1 - (0.05 / 2), n - 1) * se,
#         up_ci = mean +  (1 - (0.05 / 2), n - 1) * se) %>% 
#  ungroup()

#:#:#:#::#:PLOT#:#:#:#:#:#:
plot_comp_fit_actual_ratio<-
cbind.data.frame(model2d_rate_resp_traj_df_melt_final,
                 actual=as.numeric(unlist(data15a64_rn_ratio[,"rate_resp"]))) %>% 
    ggplot(aes(x=date))+
      geom_line(aes(y=mean), size=1, color="darkblue") +
      geom_line(aes(y=actual), size=1, color="black")+
  #geom_ribbon(aes(x=date, y=mean, ymin = lo_ci, ymax = up_ci), fill = "grey70")+
  geom_ribbon(aes(ymin=q025,ymax=q975),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  ylab("Rate of Hospitalizations per Consultations per 1,000 population")+
  xlab("Date")+
   theme(strip.text.x = element_text(size = 8, face = "bold"),
        legend.position = "bottom",
        plot.caption=element_text(hjust = 0),
        strip.background  = element_blank(),
        axis.text.x=element_text(angle = -90, hjust = 0),
        legend.title = element_blank())
  #scale_y_continuous(labels = scales::percent)
   # scale_x_date(labels = scales::date_format("%Y-%m"),
  #             limits = c(as.Date("2015-01-01"),as.Date("2019-12-31")),
  #               breaks = "3 months")
  #facet_zoom(xlim = c(253, nrow(actual_vs_fitted_and_forecast_arima)))
plot_comp_fit_actual_ratio+
    ggforce::facet_zoom(x = date>=as.Date("2019-10-01") & date<=as.Date("2019-12-31"),shrink=F,zoom.size = .5)+
    geom_vline(xintercept = as.Date("2019-10-21"), col = 2, lty = 2, size=1)+
  labs(caption="Note. Vertical Line, Social Protests;\nBlack Line= Actual Trend; Blue Line= Average Estimations")
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

vector_bsts_models<-c("impact3d1_cons_trauma", "impact3d_cons_resp","impact3d_hosp_trauma", "impact3d1_hosp_resp", "impact3d_ratio", "impact3d_ratio_resp")
names_outcomes<-c("Trauma Consultations","Respiratory Consultations","Trauma Hospitalizations","Respiratory Hospitalizations","Trauma Hospitalizations per 1,000 consultations","Respiratory Hospitalizations per 1,000 consultations")
df_results_bsts<-data.frame()
for (i in 1:6) {
  x<-vector_bsts_models[i]
dt_bsts<-cbind(
  outcome=names_outcomes[i],
  AE=sprintf("%4.2f",round(get(x)$summary$AbsEffect[1],2)),
  IC95_AE=paste0(sprintf("%4.2f",round(get(x)$summary$AbsEffect.lower[1],2)),", ",sprintf("%4.2f",round(get(x)$summary$AbsEffect.upper[1],2))),
  p=sprintf("%5.3f",round(get(x)$summary$p[1],5)),
  RE=sprintf("%4.2f",round(get(x)$summary$RelEffect[1]*100,2)),
  IC95_RE=paste0(sprintf("%4.2f",round(get(x)$summary$RelEffect.lower[1]*100,2)),", ",sprintf("%4.2f",round(get(x)$summary$RelEffect.upper[1]*100,2)))
  )
  df_results_bsts<-rbind.data.frame(df_results_bsts,dt_bsts)
}
df_results_bsts[which(df_results_bsts$p=="0.000"),"p"]<-"<0.001"

df_results_bsts %>% 
      knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption = paste0("Table 2. Estimated effects of October’s 2019 social protests on the outcomes of interest"),
               col.names = c("Outcome of Interest","Average
Effect","95%Credible Interval", "P- Value","Relative Effect(%)","95%Credible Interval"),
               align =rep('c', 101)) %>%
    kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size= 11) %>% 
  kableExtra::add_footnote(c("a- Each model had a structure of studentized distribution of errors, and a conservative prior standard deviation of .1","b- Models also included circulatory hospitalizations as a control variable.","c- Models also included circulatory consultations as a control variable.", "d- Models also included the proportion of circulatory hospitalizations of circulatory consultations (x 1,000 population) as a control variable."),
                             notation = "none")%>%
  kableExtra::scroll_box(width = "100%", height = "375px")

Table 2. Estimated effects of October’s 2019 social protests on the outcomes of interest

Outcome of Interest	Average Effect	95%Credible Interval	P- Value	Relative Effect(%)	95%Credible Interval
Trauma Consultations	-144.58	-405.43, 114.89	0.133	-15.19	-42.61, 12.08
Respiratory Consultations	-77.62	-161.06, 6.60	0.034	-45.92	-95.29, 3.90
Trauma Hospitalizations	11.95	3.39, 20.23	0.005	17.21	4.89, 29.13
Respiratory Hospitalizations	-848.49	-936.85, -756.78	<0.001	-97.64	-107.81, -87.09
Trauma Hospitalizations per 1,000 consultations	28.06	7.16, 48.16	0.006	37.93	9.68, 65.11
Respiratory Hospitalizations per 1,000 consultations	88.98	44.81, 135.78	<0.001	62.80	31.63, 95.83
a- Each model had a structure of studentized distribution of errors, and a conservative prior standard deviation of .1
b- Models also included circulatory hospitalizations as a control variable.
c- Models also included circulatory consultations as a control variable.
d- Models also included the proportion of circulatory hospitalizations of circulatory consultations (x 1,000 population) as a control variable.

2.7 Final Plots

library(cowplot)

#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-1.2 #.8
#c("impact3d1_cons_trauma", "impact3d_cons_resp","impact3d_hosp_trauma", "impact3d1_hosp_resp", "impact3d_ratio", "impact3d_ratio_resp")

Sys.setlocale(category = "LC_ALL", locale = "english")

## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"

cons_trauma_plot_red<-
plot(impact3d1_cons_trauma, "original")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=format(as.Date(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])),"%Y %b"))+
  ylim(c(0,1600))+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")

cons_resp_plot_red<-
plot(impact3d_cons_resp, "original")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=format(as.Date(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])),"%Y %b"))+
  ylim(c(0,1600))+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("b) Respiratory Consultations")

hosp_trauma_plot_red<-
plot(impact3d_hosp_trauma, "original")$data %>% 
    ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=format(as.Date(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])),"%Y %b"))+
  ylim(c(0,150))+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")

hosp_resp_plot_red<-
plot(impact3d1_hosp_resp, "original")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=format(as.Date(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])),"%Y %b"))+
  ylim(c(0,150))+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")

rate_trauma_plot_red<-
plot(impact3d_ratio, "original")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=format(as.Date(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])),"%Y %b"))+
  ylim(c(0,400))+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("e) Trauma Hospitalizations per 1,000 consultations")

rate_resp_plot_red<-
plot(impact3d_ratio_resp, "original")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=format(as.Date(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])),"%Y %b"))+
  ylim(c(0,400))+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")


#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

cons_trauma_plot<-
plot(impact3d1_cons_trauma, "original")+
    xlab("")+
    #ylab("Consultations")+
    scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])))+
    ylim(c(0,1600))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")
cons_resp_plot<-
plot(impact3d_cons_resp, "original")+
    xlab("")+
    #ylab("Consultations")+
    scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])))+
    ylim(c(0,1600))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("b) Respiratory Consultations")

hosp_trauma_plot<-
plot(impact3d_hosp_trauma, "original")+
    xlab("")+
    #ylab("Hospitalizations")+
    scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])))+
    ylim(c(0,150))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")
hosp_resp_plot<-
plot(impact3d1_hosp_resp, "original")+
    xlab("")+
    #ylab("Hospitalizations")+
    scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])))+
    ylim(c(0,150))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")

rate_trauma_plot<-
plot(impact3d_ratio, "original")+
    xlab("")+
    #ylab("Rate of Hospitalizations per Consultations")+
    ylim(c(0,400))+
    scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("e) Trauma Hospitalizations per 1,000 consultations")
rate_resp_plot<-
plot(impact3d_ratio_resp, "original")+
    xlab("")+
    #ylab("Rate of Hospitalizations per Consultations")+
    ylim(c(0,400))+
    scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")

plot1<-
plot_grid(
  cons_trauma_plot+theme(axis.text.x = element_blank()),
  cons_resp_plot+theme(axis.text.x = element_blank()),
  hosp_trauma_plot+theme(axis.text.x = element_blank()),
  hosp_resp_plot+theme(axis.text.x = element_blank()),
  rate_trauma_plot+theme(axis.text.x = element_blank()),
  rate_resp_plot,
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
  nrow = 6,  rel_heights = c(rep(1,5),1.52), align = "v"
#, rel_widths = c(1, 1)
)

#vertical
height_y<-13
height_x<-13
size_title<-14

plot1b<-
plot_grid(
  cons_trauma_plot+theme(axis.text.x = element_blank()),
  cons_resp_plot+theme(axis.text.x = element_blank()),
  hosp_trauma_plot+theme(axis.text.x = element_blank()),
  hosp_resp_plot+theme(axis.text.x = element_blank()),
  rate_trauma_plot,
  rate_resp_plot,#,
  nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,3),1.7)
)

plot1b_red<-
plot_grid(
  cons_trauma_plot_red+theme(axis.text.x = element_blank()),
  cons_resp_plot_red+theme(axis.text.x = element_blank()),
  hosp_trauma_plot_red+theme(axis.text.x = element_blank()),
  hosp_resp_plot_red+theme(axis.text.x = element_blank()),
  rate_trauma_plot_red,
  rate_resp_plot_red,#,
  nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,2),1.25)
)


ggdraw(plot1b_red)+
   theme(plot.background = element_rect(fill=NA, color = NA))

Figure 1. Estimated Trends of the Outcomes

dont_show=1
if(dont_show==0){
  pdf("_Fig1.pdf", height = 14, width = 23)
  ggdraw(plot1)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig1.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot1)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  pdf("_Fig1b.pdf", height = 19.5, width = 15.3)
  ggdraw(plot1b)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig1b.jpg", height = 19.5, width = 15.3, units = 'in', res = 600)
  ggdraw(plot1b)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}

if(dont_show==0){
  pdf("_Fig1c.pdf", height = 14, width = 23)
  ggdraw(plot1b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig1c.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot1b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
#add_sub(, "Note. The first panel shows the data and a counterfactual prediction for the post-treatment period (Blue dashed line);\nBlue area= Prediction intervals.", vpadding=grid::unit(0, "lines"), y = .55, x = 0.03, hjust = 0,size=9)
#plot2 <- cowplot::ggdraw(grid.arrange(p14, p21,p212,p32,p42,p52,p57, ncol = 2, nrow = 4)) + 
  # same plot.background should be in the theme of p1 and p2 as mentioned above
#  theme(plot.background = element_rect(fill=NA, color = NA))

#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-.8


cons_trauma_plot2_red<-
plot(impact3d1_cons_trauma, "pointwise")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")

cons_resp_plot2_red<-
plot(impact3d_cons_resp, "pointwise")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("b) Respiratory Consultations")

hosp_trauma_plot2_red<-
plot(impact3d_hosp_trauma, "pointwise")$data %>% 
    ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")

hosp_resp_plot2_red<-
plot(impact3d1_hosp_resp, "pointwise")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")

rate_trauma_plot2_red<-
plot(impact3d_ratio, "pointwise")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("e) Trauma Hospitalizations per 1,000 consultations")

rate_resp_plot2_red<-
plot(impact3d_ratio_resp, "pointwise")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")

vector_breaks_plot2b_red<- seq(238,nrow(data15a64_rn),3)
Sys.setlocale(category = "LC_ALL", locale = "english")

## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"

vector_dates_plot2b_red<-as.character(format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"))

plot2b_red<-
plot_grid(
  cons_trauma_plot2_red+scale_x_continuous(breaks=vector_breaks_plot2b_red,labels=vector_dates_plot2b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-800,300,200), limits = c(-800,300)),
  cons_resp_plot2_red+scale_x_continuous(breaks=vector_breaks_plot2b_red,labels=vector_dates_plot2b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-800,300,200), limits = c(-800,300)),
  hosp_trauma_plot2_red+scale_x_continuous(breaks=vector_breaks_plot2b_red,labels=vector_dates_plot2b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  hosp_resp_plot2_red+scale_x_continuous(breaks=vector_breaks_plot2b_red,labels=vector_dates_plot2b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  rate_trauma_plot2_red+scale_x_continuous(breaks=vector_breaks_plot2b_red,labels=vector_dates_plot2b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-400,400,200), limits = c(-400,400)),
  rate_resp_plot2_red+scale_x_continuous(breaks=vector_breaks_plot2b_red,labels=vector_dates_plot2b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-400,400,200), limits = c(-400,400)),
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
    nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,2),1.3)
)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

#vertical
height_y<-13
height_x<-13
size_title<-14

cons_trauma_plot2<-
plot(impact3d1_cons_trauma, "pointwise")+
    xlab("")+
    #ylab("Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")+
  theme(plot.title = element_text(size=size_title))
cons_resp_plot2<-
plot(impact3d_cons_resp, "pointwise")+
    xlab("")+
    #ylab("Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=12))+
  ggtitle("b) Respiratory Consultations")+
  theme(plot.title = element_text(size=size_title))
hosp_trauma_plot2<-
plot(impact3d_hosp_trauma, "pointwise")+
    xlab("")+
    #ylab("Hospitalizations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")+
  theme(plot.title = element_text(size=size_title))
hosp_resp_plot2<-
plot(impact3d1_hosp_resp, "pointwise")+
    xlab("")+
    #ylab("Hospitalizations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")+
  theme(plot.title = element_text(size=size_title))
rate_trauma_plot2<-
plot(impact3d_ratio, "pointwise")+
    xlab("")+
    #ylab("Rate of Hospitalizations per Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("e) Trauma Hospitalizations per 1,000 consultations")+
  theme(plot.title = element_text(size=size_title))
rate_resp_plot2<-
plot(impact3d_ratio_resp, "pointwise")+
    xlab("")+
    #ylab("Rate of Hospitalizations per Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")+
  theme(plot.title = element_text(size=size_title))

plot2<-
plot_grid(
  cons_trauma_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-800,300,200), limits = c(-800,300)),
  cons_resp_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-800,300,200), limits = c(-800,300)),
  hosp_trauma_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  hosp_resp_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  rate_trauma_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,400,100), limits = c(-100,400)),
  rate_resp_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,400,100), limits = c(-100,400)),
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
  nrow = 3, rel_widths = c(1, 1)
)

#horizontal
height_y<-15
height_x<-15
size_title<-16

plot2b<-
plot_grid(
  cons_trauma_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-800,300,200), limits = c(-800,300)),
  cons_resp_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-800,300,200), limits = c(-800,300)),
  hosp_trauma_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  hosp_resp_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  rate_trauma_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,400,100), limits = c(-100,400)),
  rate_resp_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,400,100), limits = c(-100,400)),
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
  nrow = 3, rel_widths = c(1, 1)
)

ggdraw(plot2b_red)+
   theme(plot.background = element_rect(fill=NA, color = NA))

Figure 2. Estimated Trends of the Outcomes

if(dont_show==0){
  pdf("___Fig2.pdf", width = 15.3, height = 19.5)
  ggdraw(plot2)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig2.jpg", width = 15.3, height = 19.5, units = 'in', res = 600)
  ggdraw(plot2)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  pdf("_Fig2b.pdf", height = 14, width = 23)
  ggdraw(plot2b)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig2b.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot2b)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  pdf("_Fig2c.pdf", height = 14, width = 23)
  ggdraw(plot2b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig2c.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot2b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}

#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-.8

hosp_trauma_plot3_red<-
plot(impact3d_hosp_trauma, "cumulative")$data %>% 
    ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")

hosp_resp_plot3_red<-
plot(impact3d1_hosp_resp, "cumulative")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")

cons_trauma_plot3_red<-
plot(impact3d1_cons_trauma, "cumulative")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")

cons_resp_plot3_red<-
plot(impact3d_cons_resp, "cumulative")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("b) Respiratory Consultations")

rate_trauma_plot3_red<-
plot(impact3d_ratio, "cumulative")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("e) Trauma Hospitalizations per 1,000 consultations")

rate_resp_plot3_red<-
plot(impact3d_ratio_resp, "cumulative")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")

vector_breaks_plot3b_red<- seq(238,nrow(data15a64_rn),3)
Sys.setlocale(category = "LC_ALL", locale = "english")

## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"

vector_dates_plot3b_red<-as.character(format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"))

plot3b_red<-
plot_grid(
  cons_trauma_plot3_red+scale_x_continuous(breaks=vector_breaks_plot3b_red,labels=vector_dates_plot3b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-4000,2000,1000), limits = c(-4000,2000)),
  cons_resp_plot3_red+scale_x_continuous(breaks=vector_breaks_plot3b_red,labels=vector_dates_plot3b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-4000,2000,1000), limits = c(-4000,2000)),
  hosp_trauma_plot3_red+scale_x_continuous(breaks=vector_breaks_plot3b_red,labels=vector_dates_plot3b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-250,250,100), limits = c(-250,250)),
  hosp_resp_plot3_red+scale_x_continuous(breaks=vector_breaks_plot3b_red,labels=vector_dates_plot3b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-250,250,100), limits = c(-250,250)),
  rate_trauma_plot3_red+scale_x_continuous(breaks=vector_breaks_plot3b_red,labels=vector_dates_plot3b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),
  rate_resp_plot3_red+scale_x_continuous(breaks=vector_breaks_plot3b_red,labels=vector_dates_plot3b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
    nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,2),1.25)
)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

#
#vertical
height_y<-13
height_x<-13
size_title<-14

#horizontal

hosp_trauma_plot3<-
plot(impact3d_hosp_trauma, "cumulative")+
    xlab("")+
   # ylab("Hospitalizations")+
      xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")+
  theme(plot.title = element_text(size=size_title))
hosp_resp_plot3<-
plot(impact3d1_hosp_resp, "cumulative")+
    xlab("")+
   # ylab("Hospitalizations")+
      xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")+
  theme(plot.title = element_text(size=size_title))
cons_trauma_plot3<-
plot(impact3d1_cons_trauma, "cumulative")+
    xlab("")+
   # ylab("Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")+
  theme(plot.title = element_text(size=size_title))
cons_resp_plot3<-
plot(impact3d_cons_resp, "cumulative")+
    xlab("")+
   # ylab("Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("b) Respiratory Consultations")+
  theme(plot.title = element_text(size=size_title))
rate_trauma_plot3<-
plot(impact3d_ratio, "cumulative")+
    xlab("")+
   # ylab("Rate of Hospitalizations per Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("e) Trauma Hospitalizations per 1,000 consultations")+
  theme(plot.title = element_text(size=size_title))
rate_resp_plot3<-
plot(impact3d_ratio_resp, "cumulative")+
    xlab("")+
   # ylab("Rate of Hospitalizations per Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")+
  theme(plot.title = element_text(size=size_title))

plot3<-
plot_grid(
  cons_trauma_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-4000,2000,1000), limits = c(-4000,2000)),
  cons_resp_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-4000,2000,1000), limits = c(-4000,2000)),
  hosp_trauma_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-250,250,100), limits = c(-250,250)), 
  hosp_resp_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-250,250,100), limits = c(-250,250)),
  rate_trauma_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),
  rate_resp_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),#,
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
  nrow = 3, rel_widths = c(1, 1)
)

height_y<-15
height_x<-15
size_title<-16

plot3b<-
plot_grid(
  cons_trauma_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-4000,2000,1000), limits = c(-4000,2000)),
  cons_resp_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-4000,2000,1000), limits = c(-4000,2000)),
  hosp_trauma_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-250,250,100), limits = c(-250,250)), 
  hosp_resp_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-250,250,100), limits = c(-250,250)),
  rate_trauma_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),
  rate_resp_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),#,
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
  nrow = 3, rel_widths = c(1, 1)
)

ggdraw(plot3b_red)+
   theme(plot.background = element_rect(fill=NA, color = NA))

Figure 3. Estimated Trends of the Outcomes

dont_show=1
if(dont_show==0){
  pdf("___Fig3.pdf", width = 15.3, height = 19.5)
  ggdraw(plot3)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig3.jpg", width = 15.3, height = 19.5, units = 'in', res = 600)
  ggdraw(plot3)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  pdf("_Fig3b.pdf", height = 14, width = 23)
  ggdraw(plot3b)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig3b.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot3b)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  pdf("_Fig3c.pdf", height = 14, width = 23)
  ggdraw(plot3b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig3c.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot3b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}

data15a64_rn_ratio_its<-
data15a64_rn_ratio %>% 
  dplyr::filter(isoweek!=53) %>% 
  dplyr::mutate(did=factor(did))

data15a64_rn_ratio_its_trauma_hosp<-
        cbind.data.frame(
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2015),"hosp_trauma"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2016),"hosp_trauma"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2017),"hosp_trauma"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2018),"hosp_trauma"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"hosp_trauma"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"did"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"isoweek"]
        )
colnames(data15a64_rn_ratio_its_trauma_hosp)<- c("ht_2015","ht_2016","ht_2017","ht_2018","ht_2019","did","isoweek")


data15a64_rn_ratio_its_resp_hosp<-
        cbind.data.frame(
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2015),"hosp_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2016),"hosp_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2017),"hosp_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2018),"hosp_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"hosp_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"did"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"isoweek"]
        )
colnames(data15a64_rn_ratio_its_resp_hosp)<- c("hr_2015","hr_2016","hr_2017","hr_2018","hr_2019","did","isoweek")

data15a64_rn_ratio_its_trauma_cons<-
        cbind.data.frame(
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2015),"cons_trauma"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2016),"cons_trauma"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2017),"cons_trauma"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2018),"cons_trauma"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"cons_trauma"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"did"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"isoweek"]
        )
colnames(data15a64_rn_ratio_its_trauma_cons)<- c("ct_2015","ct_2016","ct_2017","ct_2018","ct_2019","did","isoweek")

data15a64_rn_ratio_its_resp_cons<-
        cbind.data.frame(
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2015),"cons_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2016),"cons_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2017),"cons_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2018),"cons_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"cons_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"did"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"isoweek"]
        )
colnames(data15a64_rn_ratio_its_resp_cons)<- c("cr_2015","cr_2016","cr_2017","cr_2018","cr_2019","did","isoweek")

data15a64_rn_ratio_its_trauma_rate<-
        cbind.data.frame(
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2015),"rate"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2016),"rate"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2017),"rate"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2018),"rate"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"rate"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"did"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"isoweek"]
        )
colnames(data15a64_rn_ratio_its_trauma_rate)<- c("rt_2015","rt_2016","rt_2017","rt_2018","rt_2019","did","isoweek")

data15a64_rn_ratio_its_resp_rate<-
        cbind.data.frame(
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2015),"rate_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2016),"rate_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2017),"rate_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2018),"rate_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"rate_resp"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"did"],
        data15a64_rn_ratio_its[which(data15a64_rn_ratio_its$year==2019),"isoweek"]
        )
colnames(data15a64_rn_ratio_its_resp_rate)<- c("rr_2015","rr_2016","rr_2017","rr_2018","rr_2019","did","isoweek")

<br>

We generated several Fixed Effects Diff in Diff regressions to adjust for unobserved unit-specific (month) and confounders (circulatory) at the same time.

<br>
  
  <div style="border: 1px solid #ddd; padding: 5px; overflow-y: scroll; height:400px; overflow-x: scroll; width:100%">

####################################################
library(plm)


data_did<-cbind(
y=c("hosp_trauma", "hosp_resp", "cons_trauma","cons_resp","rate","rate_resp"),
x=c("hosp_circ", "hosp_circ", "cons_circ","cons_circ","rate_circ","rate_circ"),
colname=c("ht_","hr_","ct_","cr_","rt_","rr_")
)


rob_se_mr <-data.frame() 
for (i in  1:6) {
  assign(paste0(data_did[i,1],"_dummy_DiD_cov"), plm::plm(as.formula(paste0(data_did[i,1],"~ did+month+",data_did[i,2])), data=data15a64_rn_ratio_its, index=c("year", "isoweek"), model="within", effect= "individual"),#fixed effects ("within")
         envir = .GlobalEnv)
  #print(paste0(i,"_mr_dummy_DiD_cov"))
  #print(get(paste0(i,"_mr_dummy_DiD_cov")))
  rob_se_mr<- dplyr::bind_rows(rob_se_mr, sqrt(diag(vcovHC(get(paste0(data_did[i,1],"_dummy_DiD_cov")),  method = "arellano"))))
}
#https://stats.stackexchange.com/questions/66973/difference-between-fixed-effects-models-in-r-plm-and-stata-xtreg

#wald<- data.frame()
#for (i in 1:6) {
#  df_wald<- cbind(data_did[i,1],
##  plm::pwaldtest(get(paste0(data_did[i,1],"_dummy_DiD_cov")), vcov = function(x) vcovHC(x, type = "HC1"))$statistic,
#  plm::pwaldtest(get(paste0(data_did[i,1],"_dummy_DiD_cov")), vcov = function(x) vcovHC(x, type = "HC1"))$parameter)
#  wald<- rbind(wald,df_wald)
#}
  
#ls()[grep("*_MR_dummy_DiD_cov*",ls())]
AICs_BICs_fes<-
cbind.data.frame(
aicbic_plm(hosp_trauma_dummy_DiD_cov,"AIC"),
aicbic_plm(hosp_resp_dummy_DiD_cov,"AIC"),
aicbic_plm(cons_trauma_dummy_DiD_cov,"AIC"),
aicbic_plm(cons_resp_dummy_DiD_cov,"AIC"),
aicbic_plm(rate_dummy_DiD_cov,"AIC"),
aicbic_plm(rate_resp_dummy_DiD_cov,"AIC"),
aicbic_plm(hosp_trauma_dummy_DiD_cov,"BIC"),
aicbic_plm(hosp_resp_dummy_DiD_cov,"BIC"),
aicbic_plm(cons_trauma_dummy_DiD_cov,"BIC"),
aicbic_plm(cons_resp_dummy_DiD_cov,"BIC"),
aicbic_plm(rate_dummy_DiD_cov,"BIC"),
aicbic_plm(rate_resp_dummy_DiD_cov,"BIC")
)
#copy_names(melt(data.frame(AICs_BICs_fes)))

library(sandwich)

hcs_models<-data.frame()
for (i in 1:6) {
hcs<-t(sapply(c("HC0", "HC1", "HC2", "HC3", "HC4"), function(x) sqrt(diag(vcovHC(get(paste0(data_did[i,1],"_dummy_DiD_cov")), type = x))))) %>% data.table(keep.rownames = T) %>% slice_max(did1) %>% dplyr::select(rn)
hcs_models<-rbind(hcs_models,hcs)
}
#HC4 - small samples with influential observations HAC - heteroskedasticity and autocorrelation consistent (type ?vcovHAC for more details)
sjPlot::tab_model(
  list(hosp_trauma_dummy_DiD_cov,
                       hosp_resp_dummy_DiD_cov,
                       cons_trauma_dummy_DiD_cov,
                       cons_resp_dummy_DiD_cov,
                       rate_dummy_DiD_cov,
                       rate_resp_dummy_DiD_cov),
                  terms = c("did1"), 
                  vcov.fun = "HC",
                  string.est = "B",
                  collapse.ci = T,
                  robust= F,
                  digits= 2,
                  digits.p= 3,
                  show.aic = T,
                  ci.hyphen= ";",
                  p.style = "numeric_stars",
                  vcov.type = "HC4",
                  show.loglik= T,
                  show.p=T, 
                  show.df=F,
                  pred.labels= c("Social Protests"),
                  dv.labels= c("Trauma Hosp.","Respiratory Hosp.",
                                 "Trauma Cons.", "Respiratory Cons.",
                                 "Trauma Hosp. of Cons.","Respiratory Hosp. of Cons."), 
                  title="Estimated effect of the Social Protests, from fixed effects difference-in-difference models")

#qqnorm(residuals(hosp_resp_dummy_DiD_cov), ylab = 'Residuals')
#qqline(residuals(hosp_resp_dummy_DiD_cov))
#qqnorm(residuals(hosp_resp_dummy_DiD_cov), ylab = 'Residuals')
#qqline(residuals(hosp_resp_dummy_DiD_cov))

</div>

fixef(hosp_trauma_dummy_DiD_cov)
summary(fixef(hosp_trauma_dummy_DiD_cov))

fixefs <- merge(data15a64_rn_ratio_its, data.frame(year = names(fixef(hosp_trauma_dummy_DiD_cov)),
fixef = as.numeric(fixef(hosp_trauma_dummy_DiD_cov))),
all.x = TRUE, by = c("year"))[ , 6]

<br>

Nevertheless, must note that these coefficients are vulnerable to serial correlation. This is why difference-in-difference models were calculated in Stata, using the `xtscc` module. The code is accessible below.

<br>

Resulting data is available in the following link (https://drive.google.com/uc?export=download&id=1_wamggiWaULbFbolokBlwmr2TZOSwubN).

cap ssc install outreg2
cap install xtscc
cap ssc install coefplot
cap ssc install xttest3
cap ssc install xttest2
cap ssc install xtcsd

******************
clear all
******************
if _rc == 0 { 
    if `"`c(os)'"' == "Windows" & `"`c(username)'"' == "andre" global sf `"G:\Mi unidad\linkabbddyscriptderpaperestallidosocial"'
    cd `"${sf}"'
} 
else if _rc == 0 { 
    if `"`c(os)'"' == "Windows" & `"`c(username)'"' == "CISS Fondecyt" global sf `"G:\Mi unidad\Alvacast\CURES2_DB"'
    cd `"${sf}"'
    } 
    else if _rc == 0 { 
    if `"`c(os)'"' == "MacOSX" global sf `"/volumes/sdrive/data//"'
    cd `"${sf}"'
    }
    else if _rc == 0 { 
    mkdir "~\Stata_data" 
    } 
    else rmdir "~\Stata_data" 
 global sf `"~\Stata_data"' //if does not work
 *set working directory
 cd `"${sf}"'
 
    di `"${sf}"' 
    pwd

cap copy "https://drive.google.com/uc?export=download&id=1_wamggiWaULbFbolokBlwmr2TZOSwubN" data15a64_rn_ratio_its_did.dta,replace
use data15a64_rn_ratio_its_did.dta
*use data15a64_rn_ratio_its_did.dta

xtset year isoweek

replace did=0 if year!=2019
replace did=0 if year==2019 & isoweek<43
generate byte didf=recode(did,0,1)
drop did
gen did= didf

*semean
program define semean, rclass byable(recall) sortpreserve
    version 9.0
    syntax varlist(max=1 ts numeric) [if] [in] ///
        [, noPRInt FUNCtion(string)]
    marksample touse
    tempvar target
    if "`function'" == "" {
        local tgt "`varlist'"
    }
    else {
        local tgt "`function'(`varlist')"
    }
    capture tsset
    capture generate double `target' = `tgt' if `touse'
    if _rc > 0 {
        display as err "Error: bad function `tgt'"
        error 198
        }
    quietly summarize `target' 
    scalar semean = r(sd)/sqrt(r(N))
    if ("`print'" ~= "noprint") {
        display _n "Mean of `tgt' = " r(mean) ///
        " S.E. = " semean
    }
    return scalar semean = semean
    return scalar mean = r(mean)
    return scalar N = r(N)
    return local var `tgt'
end

*******************************************************************************


* Graficar tendencias paralelas:

local name `" "Trauma,Hospitalizations" "Respiratory,Hospitalizations" "Trauma,Consultations" "Respiratory,Consultations" "Trauma Hospitalizations,per Consultations (x1000)" "Respiratory Hospitalizations,per Consultations (x1000)" "' 
        foreach v of varlist hosp_trauma hosp_resp ///
                             cons_trauma cons_resp ///
                             rate rate_resp {
            gettoken item it : it
            gettoken nam name : name

        tabulate isoweek tx, s(`v') means
        bys isoweek: egen `v'_mean_trat = mean(`v' * (tx==1) )
        bys isoweek: egen `v'_mean_contr = mean(`v' * (tx==0) )
        g `v'_dif_tend = `v'_mean_trat - `v'_mean_contr
        label var `v'_mean_trat "`nam'"
        label var `v'_mean_contr "Average `nam'"
}

        foreach v of varlist hosp_trauma hosp_resp ///
                             cons_trauma cons_resp ///
                             rate rate_resp {
                    gettoken nam name : name                     
tw scatter `v'_mean_trat `v'_mean_contr isoweek || ///
    lfit `v'_mean_trat isoweek if tx==1 & txtime==0 || ///
    lfit `v'_mean_trat isoweek if tx==1 & txtime==1  || ///
    lfit `v'_mean_contr isoweek if tx==0 & txtime==0 || ///
    lfit `v'_mean_contr isoweek if tx==0 & txtime==1 , ///
    xlabel(1(4)52) ///
    msize(vsmall) ///
    legend(region(lstyle(none)col(none)) pos(1) ring(0) col(1) symysize(zero) /// 
    keygap(1) symxsize(large) order(1 2 3 4) lab(1 "2019") lab(2 "2015-2018") lab(3 "2019 Pre Int") lab (4 "2019 Post Int") size(tiny)) ///
    ytitle("`nam'", size(tiny)) ///
    tline(42, lpattern("_") lwidth(1) lcolor(red*0.15)) ///
    graphregion(color(gs16)) 
    graph save `v'_mean_interrupt,replace
}
graph combine hosp_trauma_mean_interrupt.gph hosp_resp_mean_interrupt.gph, graphregion(color(gs16)) 
graph combine cons_trauma_mean_interrupt.gph cons_resp_mean_interrupt.gph, graphregion(color(gs16)) 
graph combine rate_mean_interrupt.gph rate_resp_mean_interrupt.gph, graphregion(color(gs16)) 

*cons_resp_mean_interrupt cons_trauma_mean_interrupt hosp_resp_mean_interrupt ** se ven raro
*******************************************************************************
*******************************************************************************
*******************************************************************************
*COMPROBAR HETEROCEDASTICIDAD, DEPENDENCIA SERIAL, AUTOCORRELACION
*******************************************************************************
*******************************************************************************
*
local it `" "hosp_circ" "hosp_circ" "cons_circ" "cons_circ" "rate_circ" "rate_circ" "' 
local name `" "Trauma,Hospitalizations" "Respiratory,Hospitalizations" "Trauma,Consultations" "Respiratory,Consultations" "Trauma Hospitalizations,per Consultations (x1000)" "Respiratory Hospitalizations,per Consultations (x1000)" "' 
foreach v of varlist hosp_trauma hosp_resp ///
                     cons_trauma cons_resp ///
                     rate rate_resp {
                xtunitroot fisher `v', dfuller trend lags(1)
                xtunitroot fisher `v', dfuller trend lags(2)
                xtunitroot fisher `v', dfuller trend lags(3)
                xtunitroot fisher `v', dfuller trend lags(4)
                *How to do:
                ** Breusch-Pagan LM test **
                xtreg `v' i.tx##i.txtime hosp_circ i.month, fe
                xtcsd, pesaran abs

                ** Modified Wald test **
                xtreg `v' i.tx##i.txtime hosp_circ i.month, fe
                xttest3

                ** Breusch-Pagan LM test **
                xtreg `v' i.tx##i.txtime hosp_circ i.month, fe
                xttest2
}

*******************************************************************************
*******************************************************************************
*MODELO 
*******************************************************************************
*******************************************************************************

*******************************************************************************
*GENERAR WORDS MES CONTINUO Y MES FACTOR
*The Driscoll–Kraay approach provides a specific
*variant of the Newey–West robust covariance estimator computed using the Bartlett
*kernel and a time series of scores’ cross-sectional averages.
cap erase _DiD_fe.xml
local it `" "hosp_circ" "hosp_circ" "cons_circ" "cons_circ" "rate_circ" "rate_circ" "' 
local name `" "Trauma,Hospitalizations" "Respiratory,Hospitalizations" "Trauma,Consultations" "Respiratory,Consultations" "Trauma Hospitalizations,per Consultations (x1000)" "Respiratory Hospitalizations,per Consultations (x1000)" "' 
foreach v of varlist hosp_trauma hosp_resp ///
                     cons_trauma cons_resp ///
                     rate rate_resp {
    gettoken item it : it
    gettoken nam name : name
    xtscc  `v' i.did `item' month, fe 
    outreg2 using _DiD_fe, dec(2) ///
    excel e(F rss ll lag rmse p) ///
    title(Panel Estimation: Driscoll-Kraay standard errors and controlling for Circulatory Hospitalizations Consultations or Rates) ///
    stats(beta coef ci) cfmt(1) ///
    append ctitle("`nam'") nocons ///
    keep(i.did) 
}

*options ,u o ,xbu not allowed

*____________________________________________________________________________*
***GENERAR XTSCC FINAL
*____________________________________________________________________________*
*the value of consultation or hospitalizations on each value DID (including the omitted value), adjusted for the circulatory hospitalizatons/consultations/rate. 
*
*Because graphical evidence suggests
*a periodic behavior, the analysis includes the sin1 and cos1 variables, which are sine and cosine
*transformations of scaled time, respectively.

*add a fixed seasonality component based on the cosine of the season (month of year) scaled to the range (0, 2π) 
cap gen month_sin = sin(month*2*c(pi)/12)
cap gen month_cos = cos(month*2*c(pi)/12)//*(month*2*c(pi)
*gen double sintimeA= sin( 2*_pi*month/24)
*gen double costimeA= cos( 2*_pi*month/24)

**graph matrix  month_sin month_cos, ysize(4) xsize(4)
** The relation is that 2π radians equals 360◦; thus, 1 radian
**is 180◦/π or about 57.3◦. In Stata, π is wired in to as much precision as is possible in
**its calculations. T

*twoway function sine = sin(2 * pi * 12) || function cosine = cos(2 * pi * 12)
*hosp_resp ///
*cons_trauma cons_resp ///
*rate rate_resp


cap drop fe_*
cap drop xtscc_*
local it `" "hosp_circ" "hosp_circ" "cons_circ" "cons_circ" "rate_circ" "rate_circ" "' 
local name `" "Trauma,Hospitalizations" "Respiratory,Hospitalizations" "Trauma,Consultations" "Respiratory,Consultations" "Trauma Hospitalizations,per Consultations (x1000)" "Respiratory Hospitalizations,per Consultations (x1000)" "' 
foreach v of varlist hosp_trauma hosp_resp ///
                     cons_trauma cons_resp ///
                     rate rate_resp {
    gettoken item it : it
    gettoken nam name : name
    xtscc  `v' i.tx##i.txtime `item', fe
    estimates store xtscc_`v'_none
    margins, dydx(*) post 
    eststo xtscc_mar_`v'_n
    estimates restore xtscc_`v'_none
    predict fe_no_month_`v'
    xtscc  `v' i.tx##i.txtime `item' c.month, fe
    estimates store xtscc_`v'_cnt_mth
    margins i.tx##i.txtime, post
    eststo xtscc_mar_`v'_cnt_m
    estimates restore xtscc_`v'_cnt_mth
    predict fe_cont_month_`v'
    xtscc  `v' i.tx##i.txtime `item' i.month, fe
    estimates store xtscc_`v'_month
    margins i.tx##i.txtime, post
    eststo xtscc_mar_`v'_month
    estimates restore xtscc_`v'_month
    predict fe_month_`v' 
    eststo xtscc_mar_`v'_month
    xtscc  `v' i.tx##i.txtime `item' c.month#c.month, fe
    estimates store xtscc_`v'_cd_mth
    margins i.tx##i.txtime, post
    eststo xtscc_mar_`v'_cd_m
    estimates restore xtscc_`v'_cd_mth
    predict fe_cuad_`v'
    xtscc  `v' i.tx##i.txtime `item' month_cos month_sin, fe
    estimates store xtscc_`v'_sincos
    margins i.tx##i.txtime, post
    eststo xtscc_mar_`v'_sc
    estimates restore xtscc_`v'_sincos
    predict fe_sin_cos_`v'
}

*matrix list r(table) ** para ver todos los términsos
*Another way to obtain results
cap erase fe_results.csv
esttab xtscc_* ///
using fe_results.csv, append wide varlabels(1.txtime "Social Protest") keep(1.txtime) nobaselevels  ///
     stats(N N_clust r2, fmt(%9.0f %9.0f %4.3f)) ///
     cells("b(star fmt(3) label(Coef)) ci_l(fmt(2) label(CI95_Lo)) ci_u(fmt(2) label(CI95_Up)) p(fmt(%7.3f) label(p-values))") ///
     mtitles("Trauma Hospitalizations-None" "Trauma Hospitalizations-Continuous" "Trauma Hospitalizations-Factor" "Trauma Hospitalizations-Cuadratic" "Trauma Hospitalizations-SinCos" /// 
            "Respiratory Hospitalizations-None" "Respiratory Hospitalizations-Continuous" "Respiratory Hospitalizations-Factor" "Respiratory Hospitalizations-Cuadratic" "Respiratory Hospitalizations-SinCos" ///
            "Trauma Consultations-None" "Trauma Consultations-Continuous" "Trauma Consultations-Factor" "Trauma Consultations-Cuadratic" "Trauma Consultations-SinCos" ///
            "Respiratory Consultations-None" "Respiratory Consultations-Continuous" "Respiratory Consultations-Factor" "Respiratory Consultations-Cuadratic" "Respiratory Consultations-SinCos" ///
            "Trauma Hospitalizations per Trauma Consultations(x1000)-None" "Trauma Hospitalizations per Trauma Consultations(x1000)-Continuous" "Trauma Hospitalizations per Trauma Consultations(x1000)-Factor" "Trauma Hospitalizations per Trauma Consultations(x1000)-Cuadratic" "Trauma Hospitalizations per Trauma Consultations(x1000)-SinCos" ///
            "Respiratory Hospitalizations per Respiratory Consultations(x1000)-None" "Respiratory Hospitalizations per Respiratory Consultations(x1000)-Continuous" "Respiratory Hospitalizations per Respiratory Consultations(x1000)-Factor" "Respiratory Hospitalizations per Respiratory Consultations(x1000)-Cuadratic" "Respiratory Hospitalizations per Respiratory Consultations(x1000)-SinCos") ///
            legend label varwidth(25) ///
     title(Panel Estimation: Driscoll-Kraay standard errors and controlling for Circulatory Hospitalizations Consultations or Rates) ///
     compress

 *browse did hosp_trauma fe_hosp_trauma fe_month_hosp_trauma ///
        *hosp_resp fe_hosp_resp fe_month_hosp_resp ///
        *cons_trauma fe_cons_trauma fe_month_cons_trauma ///
        *cons_resp fe_cons_resp fe_month_cons_resp ///
        *rate fe_rate fe_month_rate ///
        *rate_resp fe_rate_resp fe_month_rate_resp
*cells(b(fmt(%7.4f) star label(Coef)) p(fmt(%7.3f) label(p-values))) ///
*******************************************************************************
**** CON XTREG PUEDO VER EL AIC Y EL BIC, Y LOS PUEDO COMPARAR. NO DEBIESEN SER TAN DISTINTOS CON XTSCC
*******************************************************************************

local it `" "hosp_circ" "hosp_circ" "cons_circ" "cons_circ" "rate_circ" "rate_circ" "' 
local name `" "Trauma,Hospitalizations" "Respiratory,Hospitalizations" "Trauma,Consultations" "Respiratory,Consultations" "Trauma Hospitalizations,per Consultations (x1000)" "Respiratory Hospitalizations,per Consultations (x1000)" "' 
foreach v of varlist hosp_trauma hosp_resp ///
                     cons_trauma cons_resp ///
                     rate rate_resp {
    gettoken item it : it
    gettoken nam name : name
    qui xtreg `v' i.tx##i.txtime `item', fe
    estimates store xtreg_`v'_none
    qui xtreg `v' i.tx##i.txtime `item' c.month, fe
    estimates store xtreg_`v'_cont               
    qui xtreg `v' i.tx##i.txtime `item' i.month, fe
    estimates store xtreg_`v'_fact
    qui xtreg `v' i.tx##i.txtime `item' c.month#c.month, fe
    estimates store xtreg_`v'_cuad
    qui xtreg `v' i.tx##i.txtime `item' month_sin month_cos, fe
    estimates store xtreg_`v'_sin_cos
}
*****************************************v**************************************
*****************************************v**************************************
***********TABLES************************v**************************************
qui est table xtreg_hosp_trauma_*, star b(%7.4f) stats(N r2_a aic bic rmse) keep(1.txtime)
return list
matrix stats = r(stats)'
matrix list stats
qui est table xtreg_hosp_resp_*, star b(%7.4f) stats(N r2_a aic bic rmse) keep(1.txtime)
return list
matrix stats = r(stats)'
matrix list stats
qui est table xtreg_cons_trauma_*, star b(%7.4f) stats(N r2_a aic bic rmse) keep(1.txtime)
return list
matrix stats = r(stats)'
matrix list stats
qui est table xtreg_cons_resp_*, star b(%7.4f) stats(N r2_a aic bic rmse) keep(1.txtime)
return list
matrix stats = r(stats)'
matrix list stats
qui est table xtreg_rate_none xtreg_rate_cont xtreg_rate_fact xtreg_rate_cuad  xtreg_rate_sin_cos, star b(%7.4f) stats(N r2_a aic bic rmse) keep(1.txtime)
return list
matrix stats = r(stats)'
matrix list stats
qui est table xtreg_rate_resp_*, star b(%7.4f) stats(N r2_a aic bic rmse) keep(1.txtime)
return list
matrix stats = r(stats)'
matrix list stats
*margins i.did
*nlcom (ratio1: -_b[1.did]/_b[0.did])
*****************************************


*******************************************************************************
*******************************************************************************
*******************************************************************************
*DEFINITIVE MODELS
cap erase _DiD_fe_final.xml
local it `" "hosp_circ" "hosp_circ" "cons_circ" "cons_circ" "rate_circ" "rate_circ" "' 
local term `" "month_sin month_cos" "i.month" "i.month" "i.month" "i.month" "i.month" "' 
local name `" "Trauma,Hospitalizations" "Respiratory,Hospitalizations" "Trauma,Consultations" "Respiratory,Consultations" "Trauma Hospitalizations,per Consultations (x1000)" "Respiratory Hospitalizations,per Consultations (x1000)" "' 
foreach v of varlist hosp_trauma hosp_resp ///
                     cons_trauma cons_resp ///
                     rate rate_resp {
    gettoken item it : it
    gettoken terms term : term
    gettoken nam name : name
    xtscc  `v' i.did `item' `terms', fe 
    outreg2 using _DiD_fe_final, dec(2) ///
    excel e(F rss ll lag rmse p) ///
    title(Panel Estimation: Driscoll-Kraay standard errors and controlling for Circulatory Hospitalizations Consultations or Rates) ///
    stats(beta coef ci) cfmt(1) ///
    append ctitle("`nam'") nocons ///
    keep(i.did) 
}
*******************************************************************************
*******************************************************************************



*******************************************************************************
***************************REL EFFECTS- ESP FACTOR*****************************************         
* xtscc_mar_hosp_trauma_sin_cos xtscc_mar_hosp_resp_month xtscc_mar_cons_trauma_month xtscc_mar_cons_resp_month xtscc_mar_rate_month xtscc_mar_rate_resp_month
         
estimates restore xtscc_hosp_trauma_sincos
margins i.tx##i.txtime, post
estimates replay xtscc_hosp_trauma_sincos
*# display %9.4f (65.54685-60.51798)/60.51798
*display %9.4f -3.325031/60.51798 //*LO CI
*display %9.4f 13.38276/60.51798 *UP CI
estimates restore xtscc_hosp_resp_month
margins i.tx##i.txtime, post
estimates replay xtscc_hosp_resp_month
*# display %9.4f ( 19.81461-20.44783 )/20.44783
* display %9.4f  -4.931737/20.44783
* display %9.4f  3.665299/20.44783
estimates restore xtscc_cons_trauma_month
margins i.tx##i.txtime, post
estimates replay xtscc_cons_trauma_month
*# display %9.4f (703.6545- 806.7558)/806.7558
* display %9.4f  -177.9062/806.7558
* display %9.4f -28.29627/806.7558
estimates restore xtscc_cons_resp_month
margins i.tx##i.txtime, post
estimates replay xtscc_cons_resp_month
*# display %9.4f (105.5021- 154.6745)/154.6745
* display %9.4f -79.68312/154.6745
* display %9.4f -18.66176/154.6745
estimates restore xtscc_rate_month
margins i.tx##i.txtime, post
estimates replay xtscc_rate_month
*# display %9.4f (101.063 -  76.52912)/ 76.52912
* display %9.4f 12.98471/ 76.52912
* display %9.4f  36.08302/ 76.52912
estimates restore xtscc_rate_resp_month
margins i.tx##i.txtime, post
estimates replay xtscc_rate_resp_month
*# display %9.4f (210.4664 - 137.6264)/137.6264
*  display %9.4f 38.82229/137.6264
* display %9.4f 106.8577/137.6264
*******************************************************************************
*****************************************
*****************************************


*PLOT OF LINEAR TRENDS
xtline hosp_trauma hosp_resp if year==2019, graphregion(color(gs16)) 
xtline cons_trauma cons_resp if year==2019,graphregion(color(gs16)) 
xtline rate rate_resp if year==2019,graphregion(color(gs16)) 

*SPECIFIC PLOTS
line hosp_trauma fe_no_month_hosp_trauma fe_cont_month_hosp_trauma fe_month_hosp_trauma fe_cuad_hosp_trauma fe_sin_cos_hosp_trauma isoweek if year==2019, ///
lcolor (cranberry navy brown erose black  teal) ///
legend(region(lstyle(none)col(none)) pos(4) ring(0) col(1) symysize(zero) /// 
keygap(1) symxsize(large) order(1 2 3 4 5 6) lab(1 "Actual") lab(2 "FE") lab(3 "FE +Month Cont") lab(4 "FE +Month Fact") lab(5 "FE +Month Cuad") lab(6 "FE +Month Cos Sin") size(vsmall)) ///
xtitle("Week No. (ISO 8601)", size(small)) ///
ytitle("Trauma Hospitalizations", size(medsmall)) ///
xlabel(1(4)52) ///
tline(42, lpattern("_") lwidth(1) lcolor(red*0.15)) ///
graphregion(color(gs16)) 

line hosp_resp fe_no_month_hosp_resp fe_cont_month_hosp_resp fe_month_hosp_resp fe_cuad_hosp_resp fe_sin_cos_hosp_resp isoweek if year==2019, ///
lcolor (cranberry navy brown erose black  teal) ///
legend(region(lstyle(none)col(none)) pos(2) ring(0) col(1) symysize(zero) /// 
keygap(1) symxsize(large) order(1 2 3 4 5 6) lab(1 "Actual") lab(2 "FE") lab(3 "FE +Month Cont") lab(4 "FE +Month Fact") lab(5 "FE +Month Cuad") lab(6 "FE +Month Cos Sin") size(vsmall)) ///
xtitle("Week No. (ISO 8601)", size(small)) ///
ytitle("Respiratory Hospitalizations", size(medsmall)) ///
xlabel(1(4)52) ///
tline(42, lpattern("_") lwidth(1) lcolor(red*0.15)) ///
graphregion(color(gs16)) 

line cons_trauma fe_no_month_cons_trauma fe_cont_month_cons_trauma fe_month_cons_trauma fe_cuad_cons_trauma fe_sin_cos_cons_trauma isoweek if year==2019, ///
lcolor (cranberry navy brown erose black  teal) ///
legend(region(lstyle(none)col(none)) pos(4) ring(0) col(1) symysize(zero) /// 
keygap(1) symxsize(large) order(1 2 3 4 5 6) lab(1 "Actual") lab(2 "FE") lab(3 "FE +Month Cont") lab(4 "FE +Month Fact") lab(5 "FE +Month Cuad") lab(6 "FE +Month Cos Sin") size(vsmall)) ///
xtitle("Week No. (ISO 8601)", size(small)) ///
ytitle("Trauma Consultations", size(medsmall)) ///
xlabel(1(4)52) ///
tline(42, lpattern("_") lwidth(1) lcolor(red*0.15)) ///
graphregion(color(gs16)) 

line cons_resp fe_no_month_cons_resp fe_cont_month_cons_resp fe_month_cons_resp fe_cuad_cons_resp fe_sin_cos_cons_resp isoweek if year==2019, ///
lcolor (cranberry navy brown erose black  teal) ///
legend(region(lstyle(none)col(none)) pos(2) ring(0) col(1) symysize(zero) /// 
keygap(1) symxsize(large) order(1 2 3 4 5 6) lab(1 "Actual") lab(2 "FE") lab(3 "FE +Month Cont") lab(4 "FE +Month Fact") lab(5 "FE +Month Cuad") lab(6 "FE +Month Cos Sin") size(vsmall)) ///
xtitle("Week No. (ISO 8601)", size(small)) ///
ytitle("Respiratory Consultations", size(medsmall)) ///
xlabel(1(4)52) ///
tline(42, lpattern("_") lwidth(1) lcolor(red*0.15)) ///
graphregion(color(gs16)) 

line rate fe_no_month_rate fe_cont_month_rate fe_month_rate fe_cuad_rate fe_sin_cos_rate isoweek if year==2019, ///
lcolor (cranberry navy brown erose black  teal) ///
legend(region(lstyle(none)col(none)) pos(4) ring(0) col(1) symysize(zero) /// 
keygap(1) symxsize(large) order(1 2 3 4 5 6) lab(1 "Actual") lab(2 "FE") lab(3 "FE +Month Cont") lab(4 "FE +Month Fact") lab(5 "FE +Month Cuad") lab(6 "FE +Month Cos Sin") size(vsmall)) ///
xtitle("Week No. (ISO 8601)", size(small)) ///
ytitle("Trauma Hospitalizations per Consultations (x1000)", size(small)) ///
xlabel(1(4)52) ///
tline(42, lpattern("_") lwidth(1) lcolor(red*0.15)) ///
graphregion(color(gs16)) 

line rate_resp fe_no_month_rate_resp fe_cont_month_rate_resp fe_month_rate_resp fe_cuad_rate_resp fe_sin_cos_rate_resp isoweek if year==2019, ///
lcolor (cranberry navy brown erose black  teal) ///
legend(region(lstyle(none)col(none)) pos(2) ring(0) col(1) symysize(zero) /// 
keygap(1) symxsize(large) order(1 2 3 4 5 6) lab(1 "Actual") lab(2 "FE") lab(3 "FE +Month Cont") lab(4 "FE +Month Fact") lab(5 "FE +Month Cuad") lab(6 "FE +Month Cos Sin") size(vsmall)) ///
xtitle("Week No. (ISO 8601)", size(small)) ///
ytitle("Respiratory Hospitalizations per Consultations (x1000)", size(small)) ///
xlabel(1(4)52) ///
tline(42, lpattern("_") lwidth(1) lcolor(red*0.15)) ///
graphregion(color(gs16)) 

*******************************************************************************
*******************************************************************************
**Coefficients
*est table xtscc_hosp_trauma_none, star b(%7.4f) stats(N r2 r2_a aic bic rmse) keep(1.txtime)

coefplot (xtscc_hosp_trauma_none, asequation \ xtscc_hosp_trauma_cnt_mth, asequation \ xtscc_hosp_trauma_month, asequation \ xtscc_hosp_trauma_cd_mth, asequation \ xtscc_hosp_trauma_sincos), ///
eqlabels("None" "Continuous" "Factor" "Cuadratic" "Sine Cosine", asheadings) ///
graphregion(color(gs16)) ///
 mlabsize(small) ///
 xlabel(-10(5)30,labsize(small)) ///
 ylabel(,labsize(small)) ///
 legend(off) ///
 keep(1.txtime) ///
 xline(0) ///
 xtitle("Trauma Hospitalizations", size(small))
  
 coefplot (xtscc_hosp_resp_none, asequation \ xtscc_hosp_resp_cnt_mth, asequation \ xtscc_hosp_resp_month, asequation \ xtscc_hosp_resp_cd_mth, asequation \ xtscc_hosp_resp_sincos), ///
eqlabels("None" "Continuous" "Factor" "Cuadratic" "Sine Cosine", asheadings) ///
graphregion(color(gs16)) ///
 mlabsize(small) ///
 xlabel(-15(3)15,labsize(small)) ///
 ylabel(,labsize(small)) ///
 legend(off) ///
 keep(1.txtime) ///
 xline(0) ///
 xtitle("Respiratory Hospitalizations", size(small))
 
coefplot (xtscc_cons_trauma_none, asequation \ xtscc_cons_trauma_cnt_mth, asequation \ xtscc_cons_trauma_month, asequation \ xtscc_cons_trauma_cd_mth, asequation \ xtscc_cons_trauma_sincos), ///
eqlabels("None" "Continuous" "Factor" "Cuadratic" "Sine Cosine", asheadings) ///
graphregion(color(gs16)) ///
 mlabsize(small) ///
 xlabel(-250(50)50,labsize(small)) ///
 ylabel(,labsize(small)) ///
 legend(off) ///
 keep(1.txtime) ///
 xline(0) ///
 xtitle("Trauma Consultations", size(small))
    
coefplot (xtscc_cons_resp_none, asequation \ xtscc_cons_resp_cnt_mth, asequation \ xtscc_cons_resp_month, asequation \ xtscc_cons_resp_cd_mth, asequation \ xtscc_cons_resp_sincos), ///
eqlabels("None" "Continuous" "Factor" "Cuadratic" "Sine Cosine", asheadings) ///
graphregion(color(gs16)) ///
 mlabsize(small) ///
 xlabel(-150(25)0,labsize(small)) ///
 ylabel(,labsize(small)) ///
 legend(off) ///
 keep(1.txtime) ///
 xline(0) ///
 xtitle("Respiratory Consultations", size(small))
 
 coefplot (xtscc_rate_none, asequation \ xtscc_rate_cnt_mth, asequation \ xtscc_rate_month, asequation \ xtscc_rate_cd_mth, asequation \ xtscc_rate_sincos), ///
eqlabels("None" "Continuous" "Factor" "Cuadratic" "Sine Cosine", asheadings) ///
graphregion(color(gs16)) ///
 mlabsize(small) ///
 xlabel(0(5)40,labsize(small)) ///
 ylabel(,labsize(small)) ///
 legend(off) ///
 keep(1.txtime) ///
 xline(0) ///
 xtitle("Rate of Trauma Hospitalizations per Consultations (x1,000)", size(vsmall))

coefplot (xtscc_rate_resp_none, asequation \ xtscc_rate_resp_cnt_mth, asequation \ xtscc_rate_resp_month, asequation \ xtscc_rate_resp_cd_mth, asequation \ xtscc_rate_resp_sincos), ///
eqlabels("None" "Continuous" "Factor" "Cuadratic" "Sine Cosine", asheadings) ///
graphregion(color(gs16)) ///
 mlabsize(small) ///
 xlabel(0(25)125,labsize(small)) ///
 ylabel(,labsize(small)) ///
 legend(off) ///
 keep(1.txtime) ///
 xline(0) ///
 xtitle("Rate of Respiratory Hospitalizations per Consultations(x1,000)", size(vsmall))

3 Session Info

Sys.getenv("R_LIBS_USER")

## [1] "G:/Mi unidad/paperestallido/renv/library/R-4.0/x86_64-w64-mingw32;C:/Users/CISS Fondecyt/AppData/Local/Temp/Rtmpk3lsYJ/renv-system-library"

rstudioapi::getSourceEditorContext()

## Document Context: 
## - id:        '9ECFF96C'
## - path:      'G:/Mi unidad/paperestallido/Consolidacion_BDs_FINAL.Rmd'
## - contents:  <2671 rows>
## Document Selection:
## - [2650, 20] -- [2650, 20]: ''

      tryCatch(
                     {
      save.image(paste0(getwd(),"/","Definitive_models_2021.RData"))
      #rio::export(codebook_data, "C:/Users/andre/Dropbox/Covid-19_2020/Article_SecondManuscript/LT Environmental analysis/Databases/merged_data_post_ago.dta")
                     },
                         error = function(e){
      save.image(paste0(getwd(),"/","Definitive_models_2021.RData"))
      #rio::export(codebook_data, "C:/Users/CISS Fondecyt/Dropbox/Covid-19_2020/Article_SecondManuscript/LT Environmental analysis/Databases/merged_data_post_ago.dta")
                         }
                   )
sessionInfo()

## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
##  [1] parallel  grid      splines   stats     graphics  grDevices datasets 
##  [8] utils     methods   base     
## 
## other attached packages:
##  [1] ArCo_0.3-1          cowplot_1.0.0       gridExtra_2.3      
##  [4] ISOweek_0.6-2       MSCMT_1.3.4         SCtools_0.3.1      
##  [7] future_1.17.0       doParallel_1.0.15   iterators_1.0.12   
## [10] foreach_1.5.0       imputeTS_3.1        ggforce_0.3.2      
## [13] compareGroups_4.4.5 Hmisc_4.4-0         Formula_1.2-3      
## [16] survival_3.1-12     ggsci_2.9           astsa_1.10         
## [19] vcd_1.4-7           Epi_2.40            lmtest_0.9-37      
## [22] tsModel_0.6         foreign_0.8-80      sjPlot_2.8.4       
## [25] textreg_0.1.5       CausalImpact_1.2.4  bsts_0.9.5         
## [28] xts_0.12-0          BoomSpikeSlab_1.2.3 Boom_0.9.6         
## [31] MASS_7.3-51.6       estimatr_0.22.0     rjson_0.2.20       
## [34] janitor_2.0.1       panelView_1.1.2     zoo_1.8-8          
## [37] forecast_8.12       lattice_0.20-41     DT_0.15            
## [40] data.table_1.13.0   gganimate_1.0.6     forcats_0.5.0      
## [43] stringr_1.4.0       dplyr_1.0.1         purrr_0.3.4        
## [46] readr_1.3.1         tidyr_1.1.1         tibble_3.0.3       
## [49] tidyverse_1.3.0     htmlwidgets_1.5.1   lubridate_1.7.9    
## [52] plotly_4.9.2.1      ggplot2_3.3.2       pacman_0.5.1       
## 
## loaded via a namespace (and not attached):
##   [1] tidyselect_1.1.0        lme4_1.1-23             lpSolve_5.6.15         
##   [4] munsell_0.5.0           codetools_0.2-16        effectsize_0.3.2       
##   [7] chron_2.3-55            statmod_1.4.34          withr_2.2.0            
##  [10] colorspace_1.4-1        lpSolveAPI_5.5.2.0-17.7 highr_0.8              
##  [13] NLP_0.2-0               knitr_1.29              uuid_0.1-4             
##  [16] rstudioapi_0.11         officer_0.3.13          gbRd_0.4-11            
##  [19] listenv_0.8.0           TTR_0.23-6              labeling_0.3           
##  [22] Rdpack_0.11-1           emmeans_1.4.8           slam_0.1-47            
##  [25] polyclip_1.10-0         farver_2.0.3            vctrs_0.3.2            
##  [28] generics_0.0.2          xfun_0.16               R6_2.4.1               
##  [31] assertthat_0.2.1        scales_1.1.1            nnet_7.3-14            
##  [34] gtable_0.3.0            globals_0.12.5          timeDate_3043.102      
##  [37] rlang_0.4.7             cmprsk_2.2-10           systemfonts_0.2.3      
##  [40] lazyeval_0.2.2          acepack_1.4.1           broom_0.7.0            
##  [43] checkmate_2.0.0         yaml_2.2.1              modelr_0.1.8           
##  [46] backports_1.1.7         quantmod_0.4.17         HardyWeinberg_1.6.6    
##  [49] gridtext_0.1.1          stinepack_1.4           tools_4.0.2            
##  [52] ellipsis_0.3.1          kableExtra_1.1.0        RColorBrewer_1.1-2     
##  [55] Rsolnp_1.16             Rcpp_1.0.5              plyr_1.8.6             
##  [58] base64enc_0.1-3         progress_1.2.2          prettyunits_1.1.1      
##  [61] rpart_4.1-15            fracdiff_1.5-1          haven_2.3.1            
##  [64] cluster_2.1.0           fs_1.5.0                magrittr_1.5           
##  [67] flextable_0.5.10        reprex_0.3.0            truncnorm_1.0-8        
##  [70] mvtnorm_1.1-1           sjmisc_2.8.5            hms_0.5.3              
##  [73] evaluate_0.14           xtable_1.8-4            sjstats_0.18.0         
##  [76] jpeg_0.1-8.1            etm_1.1                 readxl_1.3.1           
##  [79] shape_1.4.5             ggeffects_0.15.1        compiler_4.0.2         
##  [82] mice_3.11.0             writexl_1.3             crayon_1.3.4           
##  [85] minqa_1.2.4             htmltools_0.5.0         mgcv_1.8-31            
##  [88] ggtext_0.1.0            DBI_1.1.0               tweenr_1.0.1           
##  [91] sjlabelled_1.1.6        dbplyr_1.4.4            Rglpk_0.6-4            
##  [94] boot_1.3-25             Matrix_1.2-18           cli_2.0.2              
##  [97] quadprog_1.5-8          insight_0.9.0           pkgconfig_2.0.3        
## [100] numDeriv_2016.8-1.1     xml2_1.3.2              webshot_0.5.2          
## [103] estimability_1.3        bibtex_0.4.2.2          rvest_0.3.6            
## [106] snakecase_0.11.0        digest_0.6.25           parameters_0.8.2       
## [109] tm_0.7-7                rmarkdown_2.6           cellranger_1.1.0       
## [112] htmlTable_2.0.1         gdtools_0.2.2           curl_4.3               
## [115] urca_1.3-0              nloptr_1.2.2.2          lifecycle_0.2.0        
## [118] nlme_3.1-148            jsonlite_1.7.0          tseries_0.10-47        
## [121] viridisLite_0.3.0       fansi_0.4.1             pillar_1.4.6           
## [124] httr_1.4.2              glue_1.4.1              bayestestR_0.7.2       
## [127] zip_2.1.1               png_0.1-7               glmnet_4.0-2           
## [130] stringi_1.4.6           performance_0.4.8       blob_1.2.1             
## [133] latticeExtra_0.6-29     renv_0.12.5